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Abstract 

In order to avoid the state space explosion problem encountered in the 
quantitative analysis of large scale PEPA models, a fluid approximation ap- 
proach has recently been proposed, which results in a set of ordinary dif- 
ferential equations (ODEs) to approximate the underlying continuous time 
Markov chain (CTMC). This paper presents a mapping semantics from PEPA 
to ODEs based on a numerical representation scheme, which extends the class 
of PEPA models that can be subjected to fluid approximation. Furthermore, 
we have established the fundamental characteristics of the derived ODEs, 
such as the existence, uniqueness, boundedness and nonnegativeness of the 
solution. The convergence of the solution as time tends to infinity for several 
classes of PEPA models, has been proved under some mild conditions. For 
general PEPA models, the convergence is proved under a particular condi- 
tion, which has been revealed to relate to some famous constants of Markov 
chains such as the spectral gap and the Log-Sobolev constant. This thesis 
has established the consistency between the fluid approximation and the un- 
derlying CTMCs for PEPA, i.e. the limit of the solution is consistent with the 
equilibrium probability distribution corresponding to a family of underlying 
density dependent CTMCs. 
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1. Introduction 

Stochastic process algebras, such as PEPA [J, TIPP Q, EMPA 
IMC |H, are powerful modelling formalisms for concurrent systems which 
have enjoyed considerable success over the last decade. Such modeling can 
help designers and system managers by allowing aspects of a system which 
are not readily tested, such as scalability, to be analysed before a system 
is deployed. However, both model construction and analysis can be chal- 
lenged by the size and complexity of large scale systems. This problem, 
called state space explosion, is inherent in the discrete state approach em- 
ployed in stochastic process algebras and many other formal modelling ap- 
proaches. To overcome this problem, many work devoted exploiting the 
compositionality of theprocess algebra to decompose or simplify the under- 
lying CTMC e.g.0, i H S, B [l~o| . Another technique is the use of Kronecker 



algebra UB 3 lU JQ In addition, abstract Markov chains and stochastic 
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bounds techniques have also been used to analyse large scale PEPA mod- 



els 15, 16 



The techniques reported above are based on the discrete state space. 
Therefore as the size of the state space is extremely large, these techniques are 
not always strong enough to handle the state space explosion problem. For 
example, in the modelling of biochemical mechanisms using the stochastic 



7r-calculus (17|, |18| and PEPA (19|, |20|, the state space explosion problem 



becomes almost insurmountable. Consequently in many cases models are 
analysed by discrete event simulation rather than being able to abstractly 
consider all possible behaviours. 

To avoid this problem Hillston proposed a radically different approach 



m [2l| from the following two perspectives: choosing a more abstract state 
representation in terms of state variables, quantifying the types of behaviour 
evident in the model; and assuming that these state variables are subject 
to continuous rather than discrete change. This approach results in a set of 
ODEs, leading to the evaluation of transient, and in the limit, steady state 
measures. 

However, there are not many discussions on the fundamental problems, 
such as the existence, uniqueness, boundedness and nonnegativeness of the 
solution, as well as the its asymptotic behaviour as time tends to infinity, and 
the relationship between the derived ODEs and the underlying CTMCs for 
general PEPA models. Solving these problems can not only bring confidence 
in the new approach, but can also provide new insight into, as well as a 
profound understanding of, performance formalisms. This paper will focus 
on these topics and give answers to these problems. 

The remainder of this paper is structured as follows. Section 2 will give a 
brief introduction to PEPA as well as a numerical representation scheme de- 
veloped for PEPA. Based on this scheme, the fluid approximation of PEPE 
models will be introduced in Section 3. In this section, the existence and 
uniqueness of solutions of the derived ODEs will be presented. In addition, 
we will show that for a PEPA model without synchronisations, the solution 
of the ODEs converges as time goes to infinity and the limit coincides with 
the steady-state probability distribution of the underlying CTMC. In Sec- 
tion 4, we demonstrate the consistency between the fluid approximation and 
the Markov chains underlying the same PEPA model. This relationship will 
be utilised to investigate the long-time behaviour of the ODEs' solutions in 
Section 5. The convergence of the solutions will be proved under a partic- 
ular condition, which relates the convergence problem to some well-known 
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constants of Markov chains such as the spectral gap and the Log-Sobolev 
constant. Section 6 and 7 present an analytic approach to analyse the fluid 
approximation. For several classes of PEPA models, the convergence will be 
demonstrated under some mild conditions, and the coefficient matrices of the 
derived ODEs have been exposed to have the following property: all eigen- 
values are either zeros or have negative real parts. In addition, the structural 
property of invariance in PEPA models will be shown to play an important 
role in the proof of convergence. Finally, after presenting some related work 
in Section 8, we conclude the paper in Section 9. 



2. The PEPA modelling formalism 

This section will briefly introduce the PEPA language and its numeri- 
cal representation scheme. The numerical representation scheme for PEPA 



was developed by Ding in his thesis [22|, and represents a model numeri 



cally rather than syntactically supporting the use of mathematical tools and 
methods to analyse the model. 

2.1. Introduction to PEPA 

PEPA (Performance Evaluation Process Algebra) [lj], developed by Hill- 
ston in the 1990s, is a high-level model specification language for low-level 
stochastic models, and describes a system as an interaction of the compo- 
nents which engage in activities. In contrast to classical process algebras, 
activities are assumed to have a duration which is a random variable gov- 
erned by an exponential distribution. Thus each activity in PEPA is a pair 
(a, r) where a is the action type and r is the activity rate. The language has 
a small number of combinators, for which we provide a brief introduction be- 
low; the structured operational semantics can be found in The grammar 
is as follows: 

S ::= (a,r).S\S + S\C s 
P :;= P\X\P\P/L\C 

where S denotes a sequential component and P denotes a model component 
which executes in parallel. C stands for a constant which denotes either a 
sequential component or a model component as introduced by a definition. 
Cs stands for constants which denote sequential components. The effect of 
this syntactic separation between these types of constants is to constrain 
legal PEPA components to be cooperations of sequential processes. 
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Prefix: The prefix component (a,r).S has a designated first activity 
(a,r), which has action type a and a duration which satisfies exponential 
distribution with parameter r, and subsequently behaves as S. 

Choice: The component S + T represents a system which may behave 
either as S or as T. The activities of both S and T are enabled. Since each 
has an associated rate there is a race condition between them and the first 
to complete is selected. This gives rise to an implicit probabilistic choice 
between actions dependent of the relative values of their rates. 

Hiding: Hiding provides type abstraction, but note that the duration of 
the activity is unaffected. In P/L all activities whose action types are in L 
appear as the "private" type r. 

Cooperation:P 1X1 Q denotes cooperation between P and Q over action 
types in the cooperation set L. The cooperands are forced to synchronise on 
action types in L while they can proceed independently and concurrently with 
other enabled activities (individual activities). The rate of the synchronised 
or shared activity is determined by the slower cooperation (see [lj for details). 
We write P \\ Q as an abbreviation for P 1X1 Q when L = and P[N] is used 
to represent N copies of P in a parallel, i.e. P[3] — P \\ P || P. 

Constant: The meaning of a constant is given by a defining equation 
such as A = P. This allows infinite behaviour over finite states to be defined 
via mutually recursive definitions. 

On the basis of the operational semantic rules (please refer to [l[ for 
details), a PEPA model may be regarded as a labelled multi-transition system 



where C is the set of components, Act is the set of activities and the multi- 

relation —4 is given by the rules. If a component P behaves as Q after it 

completes activity (ot,r), then denote the transition as P^^Iq. 

The memoryless property of the exponential distribution, which is satis- 
fied by the durations of all activities, means that the stochastic process un- 
derlying the labelled transition system has the Markov property. Hence the 
underlying stochastic process is a CTMC. Note that in this representation 
the states of the system are the syntactic terms derived by the operational 
semantics. Once constructed the CTMC can be used to find steady state or 
transient probability distributions from which quantitative performance can 
be derived. 
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2.2. Numerical Representation of PEPA Models 

As explained above there have been two key steps in the use of fluid 
approximation for PEPA models: firstly, the shift to a numerical vector 
representation of the model, and secondly, the use of ordinary differential 
equations to approximate the dynamic behaviour of the underlying CTMC. 
In this paper we are only concerned with the former modification. 

This section presents the numerical representation of PEPA models de- 



(a,r) 



veloped in |22[. For convenience, we may represent a transition P —-4 Q as 



(a r P->< ^) a 

P '-^-t- Q, or often simply as P — > Q since the rate is not pertinent to 



structural analysis, where P and Q are two local derivatives. Following [21 
hereafter the term local derivative refers to the local state of a single sequen- 
tial component. In the standard structured operational semantics of PEPA 
used to derive the underlying CTMC, the state representation is syntactic, 
keeping track of the local derivative of each component in the model. In the 
alternative numerical vector representation some information is lost as states 
only record the number of instances of each local derivative: 

Definition 1. (Numerical Vector Form \2&] ). For an arbitrary PEPA 
model M. with n component types Ci,i = 1,2, ••• ,n, each with di distinct 
local derivatives, the numerical vector form of M., m(A4), is a vector with 
d = XT=i entries- The entry m[CVj records how many instances of the jth 
local derivative of component type C{ are exhibited in the current state. 

In the following we will find it useful to distinguish derivatives according 
to whether they enable an activity, or are the result of that activity: 

Definition 2. (Pre and post local derivative) 

1. If a local derivative P can enable an activity a, that is P ■, then P 
is called a pre local derivative of a. The set of all pre local derivatives 
of a is denoted by pre(a), called the pre set of a. 

2. If Q is a local derivative obtained by firing an activity a, i.e. --^Q, 
then Q is called a post local derivative of a. The set of all post local 
derivatives is denoted by post (a), called the post set of a. 

3. The set of all the local derivatives derived from P by firing a, i.e. 

post (P, a) = {Q | P^Q}, 
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is called the post set of a from P. 

Within a PEPA model there may be many instances of the same activity 
type but we will wish to identify those that have exactly the same effect 
within the model. In order to do this we additionally label activities accord- 
ing to the derivatives to which they relate, giving rise to labelled activities: 

Definition 3. (Labelled Activity). 

1. For any individual activity a, for each P G pre(a),Q G post (P, a), 
label a as a p_ ^. 



2. For a shared activity a, for each (Qi, Q2, ■ ■ • , Qk) in 

post(pre(a)[l], a) x post(pre(a)[2], a) x ••• x post(pre(a)[/c], a), 

label a as ct w , where 

w = (pre(a)[l] Q u pre(a)[2] Q 2 , ■ ■ • ,pre(a)[/c] ->■ Q k ). 

Each a p ~^® or ct w is called a labelled activity. The set of all labelled 
activities is denoted by Afa^i ■ For the above labelled activities a p_5> ^ and 
a w , their respective pre and post sets are defined as 

pre(a p ^ Q ) = {P}, post(a p ^ Q ) = {Q}, 

pre(a w ) = pre(a), post(a UI ) = {Q x , Q 2 ,--- , Qk}- 

In the numerical representation scheme, the transitions between states of 
the model are represented by a matrix, termed the activity matrix — this 
records the impact of the labelled activities on the local derivatives. 



Definition 4. (Activity Matrix, Pre Activity Matrix, Post Activity 
Matrix). For a model 1 
derivatives, the activity 1 
are defined as following 



Matrix). For a model with Nj,. ^ ^ labelled activities and Nt> distinct local 
derivatives, the activity matrix C is an Nx> x ^ ^ matrix, and the entries 



+ 1 if Pi G post(<x^ 

X p i,aj) = { -1 if Pit pre^) 
otherwise 
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where ctj is a labelled activity. The pre activity matrix C pre and post activity 
matrix C post are defined as 

C F ™(P l ,a ] ) = { + 1 ^pre(a,) > 

v ' 31 \ otherwise. ' 

QPost/p. x = f +1 if Pit P0St(aj) 

J otherwise. 

From Definitions [3] and HI each column of the activity matrix corresponds 
to a system transition and each transition can be represented by a column 
of the activity matrix. The activity matrix equals the difference between 
the pre and post activity matrices, i.e. C = C Post — C Pre . The rate of the 
transition between states is specified by a transition rate function, but we 
omit this detail here since we are concerned with qualitative analysis. See 



22] for details. 

We first give the definition of the apparent rate of an activity in a local 
derivative. 

Definition 5. (Apparent Rate of a in P) Suppose a is an activity of a 
PEP A model and P is a local derivative enabling a (i.e. P £ pre(a) y ). Let 
post (P, a) be the set of all the local derivatives derived from P by firing a, 

i.e. post(P,a) = {Q\P {a ^ Q) Q}. Let 

r a (P)= £ r£-* (1) 

QSpost(P,a) 

The apparent rate of a in P in state x ; denoted by r a (x, P), is defined as 

r a (x,P)=x[P]r a (P). (2) 

The above definition is used to define the following transition rate func- 
tion. 

Definition 6. (Transition Rate Function) Suppose a is an activity of a 
PEPA model and x denotes a state vector. 

1. If a is individual, then for each P — > Q, the transition rate func- 
tion of labelled activity a p ^ Q in state x is defined as 

/(x,a p ^)=x[P]r^. (3) 
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2. If a is synchronised, with pre(a) = {Pi, P 2 , • • • , Pk], then for each 

(Qi, Q2, ■ ■ ■ , Qk) e post(Pi, a) x post(P 2 , a) x • • • x post(P fc , a), 

let w = (Pi — > Qi, P2 — > Q2, ■ ■ • , Pfc — > Qk)- Then the transition rate 
function of labelled activity a w in state x is defined as 

( * r Pi->Qi \ 

\~1 r a(P) J ie{l,-,k} 

where r a (x, Pj) = x[Pi]r a (Pj) is the apparent rate of a in Pi in state x. 
So 

Note that Definition [6] accommodates the passive or unspecified rate T. 
An algorithm for automatically deriving the numerical representation of a 



PEPA model was presented in [22 



Remark 1. Definition^ accommodates the passive or unspecified rate T. If 
there are some rj 7-5 ^ which are T , then the relevant calculation in the rate 
functions and @) can be made according to the following inequalities and 
equations that define the comparison and manipulation of unspecified activity 
rates (see Section 3.3.5 in fij): 



r < wT for all r G R and for all w G N 

w±T < w 2 T if w\ < W2 for all W\, w 2 G N 

w{T + w 2 T = (wi + W2)T for all W\, w 2 G N 

Moreover, we assume that ■ T = 0. So the terms such as 'min{AT, rB} " 
are interpreted as f^ij: 



minjylT, rB} 



rB, A > 0, 
0, A = 0. 



The transition rate function has the following properties [22 
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Proposition 1. The transition rate function is nonnegative; if P is a pre 
local derivative of I, i.e. P G pre(Z), then the transition rate function of I in 
a state x is less than the apparent rate of I in U in this state, that is 



o</(x,z)<n(x,P) = x[P] n (P), 

where ri(P) is the apparent rate of I in P for a single instance of P. 

Proposition 2. Let I be an labelled activity, and x, y be two states. The 
transition rate function /(x, I) defined in Definition^ satisfies: 

1. For any H > 0, Hf(x/H,l) = /(x,f)- 

2. There exists M > such that |/(x, /) — f(y,l)\ < M||x — y|| for any 
x, y and I. 

Hereafter || ■ || denotes any matrix norm since all finite matrix norms are 
equivalent. The first term of this proposition illustrates a homogenous prop- 
erty of the rate function, while the second indicates the Lipschtiz continuous 
property, both with respect to states. 



3. Fluid approximations for PEPA models 

The section will introduce the fluid-flow approximations for PEPA mod- 
els, which leads to some kind of nonlinear ODEs. The existence and unique- 
ness of the solutions of the ODEs will be established. Moreover, a conserva- 
tion law satisfied by the ODEs will be shown. 

3.1. State space explosion problem: an illustration by a tiny example 
Let us first consider the following tiny example. 

Model 1. 

Useri =(task\, a).User 2 
User 2 =(task2,b).Useri 
Provideri =(taski, a).Provider 2 
Provider 2 =(reset,d). Provideri 
Useri\\ ■ ■ ■ \\Useri ^ Provider^- ■ -\\Provideri 
M copies N copies 
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According to the semantics of PEPA originally defined in [lj, the size of 
the state space of the CTMC underlying Model [T] is 2 M+N . That is, the size 
of the state space increases exponentially with the numbers of the users and 
providers in the system. Consequently, the dimension of the infinitesimal 
generator of the CTMC is 2 M+N x 2 M+N . The computational complexity of 
solving the global balance equation to get the steady-state probability dis- 
tribution and thus derive the system performance, is therefore exponentially 
increasing with the numbers of the components. When M and/or N are 
large, the calculation of the stationary probability distribution will be infea- 
sible due to limited resources of memory and time. The problem encountered 
here is the so-called state-space explosion problem. 

A model aggregation technique, i.e. representing the states by numerical 



vector forms, which is introduced by Gilmore et al. in [24[ and by Hillston 
in 



21], can help to relieve the state-space explosion problem. It has been 
proved in 22] that by employing numerical vector forms the size of the state 
space can be reduced to 2 M x 2 N to (M + 1) x (N + 1), without relevant 
information and accuracy loss. However, this does not imply there is no 
complexity problem. The following table, Table [TJ gives the runtimes of 
deriving the state space in several different scenarios. All experiments were 
carried out using the PEPA Plug-in (vO.0.19) for Eclipse Platform (v3.4.2), 
on a 2.66GHz Xeon CPU with 4Gb RAM running Scientific Linux 5. The 
runtimes here are elapsed times reported by the Eclipse platform. 



Tabic 1: Elapsed time of state pace derivation 



(M,N) 


(300,300) 


(350,300) 


(400,300) 


(400,400) 


time 


2879 ms 


4236 ms 


"Java heap space" 


"GC overhead limit 
exceeded" 



If there are 400 users and 300 providers in the system, the Eclipse platform 
reports the error message of "Java heap space", while 400 users and 400 
providers result in the error information of "GC overhead limit exceeded". 
These experiments show that the state-space explosion problem cannot be 
completely solved by just using the technique of numerical vector form, even 
for a tiny PEPA model. That is, in order to do practical analysis for large 
scale PEPA models we need new approaches. 



12 



3.2. Fluid approximation of PEPA models 

In the numerical representation of PEPA presented in Section 12.21 a nu- 
merical vector form is introduced to capture the state information of models 
with repeated components. In this vector form there is one entry for each 
local derivative of each component type in the model. The entries in the 
vector are no longer syntactic terms representing the local derivative of the 
sequential component, but the number of components currently exhibiting 
this local derivative. Each numerical vector represents a single state of the 
system. The rates of the transitions between states are specified by the tran- 
sition rate functions. For example, the transition from state x to x + I can 
be written as 

x — y x + /, 

where I is a transition vector corresponding to the labelled activity / (for 
convenience, hereafter each pair of transition vectors and corresponding la- 
belled activities shares the same notation), and /(x, I) is the transition rate 
function, reflecting the intensity of the transition from x to x + /. 

The state space is inherently discrete with the entries within the numerical 
vector form always being non-negative integers and always being incremented 
or decremented in steps of one. As pointed out in 21], when the numbers of 
components are large these steps are relatively small and we can approximate 
the behaviour by considering the movement between states to be continuous, 
rather than occurring in discontinuous jumps. In fact, let us consider the 
evolution of the numerical state vector. Denote the state at time t by x(t). 
In a short time At, the change to the vector x(£) will be 

x(-,t + At)-x(- ) t) = F(x(-,t))At = At U{x{;t),l). 

Dividing by At and taking the limit, At — > 0, we obtain a set of ordinary 
differential equations (ODEs): 

I - *M, <*) 

i/(x,(). (6) 



where 



^£-4label 

Once the activity matrix and the transition rate functions are generated, 
the ODEs are immediately available. All of them can be obtained automat- 



ically by a derivation algorithm presented in [22 
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Let U be a local derivative. For any transition vector /, l[U] is either ±1 
or 0. If l[U] = — 1 then U is in the pre set of /, i.e. U G pre(Z), while l[U] = 1 
implies U G post(Z). According to (J5J) and (jSJ), 



dx(Z7, t) v^ 7rrn ,/ n 
— 3} — = 2^ l W\f{^ 1 ) 



d* 

= - £ /(x,/)+ £ /(x,/) (?) 

W[£/]=-l M[£/]=l 

= - £ /(x,0+ £ /(x,0- 
0|[/eprc(Z)} {i|c/e P ost(/)} 

The term X^{;|c/epre(0} ^( x ' ^ represents the "exit rates" in the local derivative 
U, while the term X]{;|c/epost(«)} ^( x ' re fl ec ts the "entry rates" in U. The 
formulae (jSJ) and ([6]) are activity centric while (I7j) is local derivative centric. 
Our approach to derive ODEs has extended previous results presented in 



the literature |21|, |23|, [25(, by relaxing restrictions such as allowing shared 
activities may have different local rates, each action name may appear in 
different local derivatives within the definition of a sequential component, 
and may occur multiple times with that derivative definition, etc. 

For an arbitrary CTMC, the evolution of probabilities distributed on each 
state can be described by a set of linear ODEs ([ifl, page 52). For example, 
for the (aggregated) CTMC underlying a PEPA model, the corresponding 
differential equations describing the evolution of the probability distributions 
are 

where each entry of 7r(£) represents the probability of the system being in each 
state at time t, and Q is an infinitesimal generator matrix corresponding to 
the CTMC. Clearly, the dimension of the coefficient matrix Q is the square 
of the size of the state space, which increases as the number of components 
increases. 

The derived ODEs (j3J) describe the evolution of the population of the 
components in each local derivative, while © reflects the the probability 
evolution at each state. Since the scale of (EJ), i.e. the number of the ODEs, 
is only determined by the number of local derivatives and is unaffected by 
the size of the state space, so it avoids the state-space explosion problem. 
But the scale of flS} depends on the size of the state space, so it suffers from 
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the explosion problem. The price paid is that the ODEs §5§ are generally 
nonlinear due to synchronisations, while flS} is linear. However, if there is 
no synchronisation contained then (jSJ) becomes linear, and there is some 
correspondence and consistency between these two different types of ODEs, 
which will be demonstrated in Section 13.41 

3. 3. Existence and uniqueness of ODEs ' solution 

For any set of ODEs, it is important to consider if the equations have a 
solution, and if so whether that solution is unique. 

Theorem 1. For a given PEP A model without passive rates, the derived 
ODEs from this model have a unique solution in the time interval [0, oo). 

Proof. Notice that each entry of F(x) = Y2i lf(x, I) is a linear combination of 
the transition rate functions /(x, /), so F(x) is globally Lipschitz continuous 
since each /(x, /) is globally Lipschitz continuous. That is, there exits M > 
such that Vx, y, 

||F(x)-F(y)||<M||x-y||. (9) 



By the classical theory in ODEs (e.g. Theorem 6.2.3 in 27j, page 14), the 



derived ODEs have a unique solution in [0, oo). □ 

As we have mentioned, in the formula (0), the term J2{i\uepvc(i)} /( x '0 
represents the exit rates in the local derivative U, while the term 
J2{i\uepost(i)} /( x ' reflects the entry rates in U. For each type of component 
at any time, the sum of all exit activity rates must be equal to the sum of 
all entry activity rates, since the system is closed and there is no exchange 
with the environment. This leads to the following proposition. 

Proposition 3. Let Ci } be a local derivative of component type Ci. Then for 

any i and t, ^ — 7 - = 0, and £\ x (C ij: t) = £\ x {C ip 0) . 

3 

Proof. In the definition of activity matrix, the numbers of —1 and 1 appearing 
in the entries of any transition vector / (i.e. a column of the activity matrix), 



which correspond to the component type Cj, are the same [22 



i.e. 



#{j : l[a) = -1} = #{j : l[a) = 1}. (10) 
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i, if/[cy = ±i, 

0, otherwise. 



y T^ =y Tj- z/(Xi0= J- y r z/(x>0 = . 
i i 



Let y be an indicator vector with the same dimension as / satisfying: 

y[CY 

So y T l = by (JTUJ) . Thus 

dx 
dt 

Edx(CV,tJ T dx / v 
— = y — = 0. So j x \pi 3 it) is a constant and equal 

j 

to 2~2j x (Cjj)O)) i- e - the number of the copies of component type in the 
system initially. □ 

Proposition [3] means that the ODEs satisfy a Conservation Law, i.e. the 
number of each kind of component remains constant at all times. 

3-4- Convergence and consistence of ODEs' solution: nonsynchronised mod- 
els 

Now we consider PEPA models without synchronisation. For this spe- 
cial class of PEPA models, we will show that the solutions of the derived 
ODEs have finite limits. Moreover, the limits coincide with the steady-state 
probability distributions of the underlying CTMCs. 

3.4-1- Features of ODEs without synchronisations 

Suppose the PEPA model has no synchronisations. Without loss of gener- 
ality, we suppose that there is only one kind of component C in the system. 
In fact, if there are several types of component in the system, the ODEs 
related to the different types of component can be separated and treated in- 
dependently since there are no interactions between them. Thus, we assume 
there is only one kind of component C and that C has k local derivatives: 
C X ,C 2 ,>- ,C k . Then (ED is 

d(x(Ci,t),^,x(C fc ,t)) r = ^^ (x> ^ (n) 
dt i 
Since (jlip are linear ODEs, we may rewrite (jlip as the following matrix form: 

d(x(C 1 ,t),---,x(C fc ,t)) T 



dt 



Q J (x(d,t),---,x(C fc ,t)r , (12) 
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where Q = (g^) is a k x k matrix. 
Q has many good properties. 

Proposition 4. Q = (qij) kxk in IfW) is an infinitesimal generator matrix, 
that is, (qij) kxk satisfies 

1. < — qu < oo for all i; 

2. qij > for all i ^ j ; 

3 - J2j=iQij = f or al1 i- 
Proof. According to (|T2|) . we have 

dX(a,t) =Ex(C,,%, (13) 



dt 



Notice by (13), 



dt 



E /( x ' z )+ E /( x '0- 

{z|Cieprc(i)} {il^epostCO} 



So 



A.- 



E x (^% = - E ^ x ' z )+ E /( x >o- ( 14 ) 

J=l 0|Ci6pre(0} 0|Ci6post(Z)} 

Since there is no synchronisation in the system, the transition func- 
tion /(x, I) is linear with respect to x and therefore there is no nonlinear 
term,"min", in it. In particular, if G pre(Z), then /(x, I) = rj(Cj)x[Ci], 
which is the apparent rate of I in d in state x defined in Definition [5j 
We should point out that according to our semantics of mapping PEPA 
models to ODEs, the fluid approximation-version of /(x, /) also holds, i.e. 
f(x(t),l) = ri(Ci)x(C h t). So (fUD becomes 

x(C i ,t)q ii + Y,x(C j ,t)q ji = x(C i ,t) E E /( x >0- 

j¥» {i|Ciepre(Z)} {Z|Ciepost(Z)} 

(15) 

Moreover, as long as /(x, /) = n(Cj)x(Cj) for some / and some positive 
constants r/(Cj), which implies that Z can be fired at Cj, we must have Cj G 
pre(/). That is to say, if Cj G post(Z) then /(x, /) cannot be of the form 
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of rx(Ci,t) for any constant r > 0. Otherwise, we have Cj G pre(Z), which 
results a contradiction^ to Cj G post(Z). So according to (JT51) . we have 

x(ci,t) fe = x(a,*) E (-n(Ci)), (16) 

0|C«6pre(0} 

JXq,, E /(x,Z). (17) 

jVi {J|Ci£post(Z)} 

Thus by (JTB]), qu = J2{i\Ci£ P re(i)} (~ r i( C i))^ and < < 00 for a11 i- Item 
1 is proved. 

Similarly, for any Cj, j ^ i, if /(x, /) = rx(Cj,t) for some / and positive 
constant r, then clearly Cj is in the pre set of I. That is Cj G pre(Z). So by 

(HID, 

{J|Cj6pre(0,Ci6post(0} * 

which implies gjj = * — f° r a ^ * 7^ ii i- e - it em 2 holds. 

We now prove item 3. By Proposition [3J, 

x(Ci,t) | x(C 2 ,t) | | x(C fc ,t) Q 
lit (it (it 

Then by (TTS]) and (USD, for all t, 

x(Ci, t) E g y + x(C 2 , t) E g 2 j H h x(C fe , t) E Qkj 

j=i j=i j=i 

= e x(c, t) e = E E x (^> ^ = E ^r 1 = °- 
i=i j=i j=i i=i j=i 

This implies X^Li %' = f° r & h □ 



1 In this paper we do not allow a self-loop in the considered model. That is, any PEPA 
finition like ' 
is not allowed. 



definition like "C = (a, r).C" which results in C € pre(a) and C € post(a) simultaneously, 
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In the proof of Proposition HJ we have shown the relationship between 
the coefficient matrix Q and the activity rates: 

qu = - Yl r i( c i)i % = ^2 r f i ^ Cj (* j)- 

{l\d£prc(l)} I 

We point out that this infinitesimal generator matrix Qkxk may not be the 
infinitesimal generator matrix of the CTMC derived via the usual semantics 
of PEPA (we call it the "original" CTMC for convenience). In fact, the 
original CTMC has a state space with k N states and the dimension of its 
infinitesimal generator matrix is k N x k , where N is the total number of 
components in the system. However, this Qkxk is the infinitesimal generator 
matrix of a CTMC underlying the PEPA model in which there is only one 
copy of the component, i.e. N — 1. To distinguish this from the original one, 
we refer to this CTMC as the "singleton" CTMC. 

3.4-2. Convergence and consistency for the ODEs 

Proposition H] illustrates that the coefficient matrix of the derived ODEs is 
an infinitesimal generator. If there is only one component in the system, then 
equation fTT2|) captures the probability distribution evolution equations of the 
original CTMC. Based on this proposition, we can furthermore determine 
the convergence of the solutions. 

Theorem 2. Suppose x(Cj,i) (j = 1,2, ■■• ,k) satisfy < f77)) . then for any 
given initial values x(Cj,0) > (j = 1,2, ■•■ ,k), there exist constants 
x(Cj,oo), such that 

]imx(C,-,f) =x(C i ,oo), j = 1,2,- •• ,jfc. (20) 

t— >oo 

Proof. By Proposition HJ the matrix Q in ( fl2|) is an infinitesimal genera- 
tor matrix. Consider a "singleton" CTMC which has the state space S = 
{Ci, C*2, • • • , Cfc}, the infinitesimal generator matrix Q in (I12p and the initial 
probability distribution vr(Cj, 0) = x< - < ^'°^ (j = 1, 2, ■ • • , k). Then according 
to Markov theory (j26|, page 52), n(Cj,t) (j = 1,2, ■■• ,k), the probability 
distribution of this new CTMC at time t, satisfies 

lW^kp^» =WCl , (),...,. (C t , f ))Q (21) 
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Since the singleton CTMC is assumed irreducible and positive-recurrent, it 
has a steady-state probability distribution {ir(Cj, oo)}^ =1 , and 

l3mn(C j ,t)=ir{C j ,oo), j = 1, 2, • • • , k. (22) 

t— >oo 

Note that x ^'^ also satisfies (12T|) with the initial values x ^'°- > equal to 
7r(Cj,0), where N is the population of the components. By the uniqueness 
of the solutions of ( l2Tj) . we have 

^M = 7r(C^), j = 1,2, (23) 

and hence by (122]) . 



lim x(Q, t) = lim Nir(Cj,t) = Nir(Cj, oo), j = 1, 2, • • ■ , k. 

t— >oo t— >oo 



□ 



Clearly, if there are multiple types of component in the system, then 
Theorem [2] holds for each each component type, since there is no cooperation 
between different component types and they can be treated independently 



It is shown in [28] that for some special examples the equilibrium solu- 
tions of the ODEs coincide with the steady state probability distributions 
of the underlying original CTMC. This theorem states that this holds for 
all for PEPA models without synchronisations. Moreover, it is also exposed 
by this theorem that the fluid approximation is consistent with the CTMC 
underlying the same nonsynchronised PEPA model. 

4. Relating to density dependent CTMCs 

A general PEPA model may have synchronisations, which result in the 
nonlinearity of the derived ODEs. However, it is difficult to rely on pure 
analytical methods to explore the asymptotic behaviour of the solution of 
the derived ODEs from an arbitrary PEPA model (except for some special 
classes of models, see the next two sections). 



Fortunately, Kurtz's theorem |29l. |30| establishes the relationship between 



a sequence of Markov chains and a corresponding set of ODEs: the complete 
solution of some ODEs is the limit of a sequence of Markov chains. In the 
context of PEPA, the derived ODEs can be considered as the limit of pure 
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jump Markov processes, as first exposed in [25( for a special case. Thus 
we may investigate the convergence of the ODEs' solutions by alternatively 
studying the corresponding property of the Markov chains through this con- 
sistency relationship. This approach leads to the result presented in the next 
section: under a particular condition the solution will converge and the limit 
is consistent with the limit steady-state probability distribution of a family 
of CTMCs underlying the given PEPA model. Let us first introduce the 
concept of density dependent Markov chains underlying PEPA models. 

4-1. Density dependent Markov chains underlying PEPA models 

In the numerical state vector representation scheme, each vector is a single 
state and the rates of the transitions between states are specified by the rate 
functions. For example, the transition from state x to x + I can be written 

as 

x —4 x + /. 

Since all the transitions are only determined by the current state rather than 
the previous ones, given any starting state a CTMC can be obtained. More 
specifically, the state space of the CTMC is the set of all reachable numerical 
state vectors x. The infinitesimal generator is determined by the transition 
rate function, 

g x ,x+i = /(x,Z). (24) 

Because the transition rate function is defined according to the seman- 
tics of PEPA, the CTMC mentioned above is in fact the aggregated CTMC 
underlying the given PEPA model. In other words, the transition rate of the 
aggregated CTMC is specified by the transition rate function in Definition [6j 

It is obvious that the aggregated CTMC depends on the starting state of 
the given PEPA model. By altering the population of components presented 
in the model, which can be done by varying the initial states, we may get a 
sequence of aggregated CTMCs. Moreover, Proposition [2] indicates that the 
transition rate function has the homogenous property: Hf(x/H, I) = /(x, I), 
WH > 0. This property identifies the aggregated CTMC to be density depen- 
dent. 

Definition 7. A family of CTMCs {X n } n is called density dependent 

if and only if there exists a continuous function /(x, /), x G R d , I G Z d ; such 
that the infinitesimal generators of X n are given by: 
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where (^x+z denotes an entry of the infinitesimal generator of X n , x a nu- 
merical state vector and I a transition vector. 

This allows us to conclude the following proposition. 

Proposition 5. Let {X n } be a sequence of aggregated CTMCs generated 
from a given PEPA model (by scaling the initial state), then {X n } is density 
dependent. 

Proof. For any n, the transition between states is determined by 

?£Ui = /(x,Z), 

where x, x + I are state vectors, / corresponds to an activity, /(x, I) is the 
rate of the transition from state x to x + /. By Proposition [21 

nf(x/n,l) = /(x,Z). 

So the infinitesimal generator of X n is given by: 

Therefore, {X n } is a sequence of density dependent CTMCs. □ 

In particular, the family of density dependent CTMCs, {X n (t)}, derived 
from a given PEPA model with the starting condition X n (0) = nxo (Vn), is 
called the density dependent CTMCs associated with x$. The CTMCs Xn ^ 
are called the concentrated density dependent CTMCs. Here n is called the 
concentration level, indicating that the entries within the numerical vector 
states (of are incremented and decremented in steps of -. 

For example, Consider the following PEPA model, which is Model [T] pre- 
sented previously: 

User 1 =(taski,a).User 2 

User 2 =(task2,b).Useri 
Provider \ ={task\, a). Provider 2 
Provider^ = (reset,d).Provideri 
[User^M]) IXI (Provider^N]) . 
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The activity matrix and transition rate functions have been specified in Ta- 
ble [2J In this table, £/j,Pj (i = 1,2) are the local derivatives representing 
Usevi and Provider \ respectively. For convenience, the labelled activities 
or transition vectors tas/ci^ c/l_5 ' c/2 ' Pl_5 ' P2 ' ) , task 2 2 ~* Ul , reset P2 ^ Pl will subse- 
quently be denoted by l taskl ; i task ? preset reS p ec tively. 



Table 2: Activity matrix and transition rate function of Model [T] 



/ 




task u 2 M 


reset P2 ^ Pl 




-1 


1 





u 2 


1 


-1 





Pi 


-1 





1 


p 2 


1 





-1 


/(x,/) 


amin(x[L r i], x[Px]) 


bx[U 2 ] 


dx[P 2 ] 



Suppose xi = (M,0,iV,0f = (1,0, 1,0) T . Let X^t) be the aggregated 
CTMC underlying Model [1] with initial state xi. Then the state space of 
Xi(t), denoted by Si, is composed of 

Xi = (1,0,1, Of, x 2 = (0,1,0, If, 

x 3 = (1,0,0, if, x 4 = (0,1,1, Of. ^ 

According to the transition rate functions presented in Table [2j we have, for 
instance, 

<?Sx 2 = Wj+j*"* = /(xi,/ iasfcl ) = amin(x 1 [t/ 1 ],x 1 [P 1 ]) = a. 

Varying the initial states we may get other aggregated CTMCs. For example, 
let X 2 {t) be the aggregated CTMC corresponding to the initial state X 2 (0) = 
2x = (2, 0, 2, Of. Then the state space S 2 of X 2 (t) has the states 

x 1 = (2,0,2,0f, x 2 = (1,1,1, If, x 3 = (l,l,2,0f, 

x 4 = (l,l,0,2f, x 5 = (0,2,1, If, x 6 = (2,0,l,lf, (26) 

x 7 = (0,2,0,2f, x 8 = (0,2,2,0) T , x 9 = (2, 0, 0, 2f . 

The rate of transition from x x to x 2 is determined by 

= ?x 1 ,x 1+ ^ 1 = f(*iJ taskl ) = 2a = 2/( Xl /2,^). 
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Similarly, let X n (t) be the aggregated CTMC corresponding to the initial 
state X n (0) = nx . Then the transition from x to x + I is determined by 

= /(x,Q = «/(*/«. 0- 

Thus a family of aggregated CTMCs, i.e. {X n (i}}, has been obtained from 
Model [TJ These derived {X n (t)} are density dependent CTMCs associated 
with x . As illustrated by this example, the density dependent CTMCs are 
obtained by scaling the starting state x . So the starting state of each CTMC 
is different, because X n (0) = nx , i.e. X n (0) = n(M, 0, N, 0) T . 

4-2. Fluid approximation as the limit of the CTMCs 

As discussed above, a set of ODEs and a sequence of density dependent 
Markov chains can be derived from the same PEPA model. The former one 
is deterministic while the latter is stochastic. However, both of them are 
determined by the same activity matrix and the same rate functions that are 
uniquely generated from the given PEPA model. Therefore, it is natural to 
believe that there is some kind of consistency between them. 

As we have mentioned, the complete solution of some ODEs can be the 



limit of a sequence of Markov chains according to Kurtz's theorem [29|, |30 



Such consistency in the context of PEPA has been previously illustrated for 



a particular PEPA model [25] . Here we give a modified version of this result 



for general PEPA models, in which the convergence is in the sense of almost 



surely rather than probabilistically as in [25 



Theorem 3. Let X(t) be the solution of the ODEs |3P derived from a given 
PEPA model with initial condition xo, and let {X n (t)} be the density de- 
pendent CTMCs associated with x underlying the same PEPA model. Let 
X n (t) = then for any t > 0, 

lim sup \\X n (u) - X(u)\\ = a.s. (27) 

u<t 



Proof. According to Kurtz's theorem, which is listed in Appendix A it is 
sufficient to prove: for any compact set K C M. Nd , 

1. 3M K > such that ||F(x) - F(y)|| < M K ||x-y||; 
2- £JKI|sup xe tf/(x,0 < oo. 
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Clearly, above term 1 is satisfied. Since /(x, /) is continuous by Proposition [21 
it is bounded on any compact K. Notice that any entry of / takes values in 
{1,-1,0}, so || 2 1| is bounded. Thus term 2 is satisfied, which completes the 
proof. □ 

Theorem [3] allows us to investigate the properties of X(t) through study- 
ing the characteristics of the family of CTMCs X n (t) = nK ' . Notice that 
X n (t) takes values in the state space which corresponds to the starting state 
nxo- Clearly, each state of X n (t) is bounded and nonnegative (more pre- 
cisely, each entry in any numerical state vector is nonnegative, and bounded 
by ||x || according to the conservative law). So the ODEs' solution X(t) 
inherits these characteristics since X(t) is the limit of X n (t) as n goes to 
infinity. That is, X(t) is bounded and nonnegative. The proof is trivial and 
omitted here. Instead, a purely analytic proof of these properties will be 
given in Section 16.11 

4-3. Consistency between the fluid approximation and the CTMCs 

As shown in Theorem |3l for a given PEPA model with synchronisations, 
the derived ODEs can be taken as the underlying density dependent CTMC 
with the concentration level infinity. If the given model has no synchroni- 
sations, then by Proposition H] and the proof of Theorem [2J the ODEs coin- 
cide the probability distribution evolution equations of the CTMC with the 
concentration level one, except for a scaling factor. These two conclusions 
embody the consistency between the fluid approximation and the CTMCs 
for PEPA models. Moreover, as we will see in the next section, the limit of 
the ODEs' solution as time tends to infinity (if it exists) is consistent with 
the limit of the expectations of the corresponding CTMCs. 

However, it is natural to ask such a question: for a PEPA model with 
synchronisations, since the ODEs corresponds the CTMC with the concen- 
tration level infinity and our performance evaluation is based on the usual 
underlying CTMC, i.e. the CTMC with the concentration level one, what is 
the loss brought by using the fluid approximation approach? Although in the 
sense of concentration level there is a gap between one and infinity, but in 
practice for large scale synchronised models (i.e. models with a large number 
of repetitive components), the relative errors of performance measure such as 
throughput and average response time between the concentration level one 



and infinity are usually small enough to be ignored (i.e. less than 5%) 22 



Therefore, it is safe to employ the fluid approximation for large scale models. 
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This paper more likely concentrates on the long-time behaviour the ODEs' 
solution. In the following sections, we focus on the problem of whether the 
ODEs' solution converges as time goes to problem will be presented in the 
next section, which is based on the consistency relationship between the 
derived ODEs and the CTMCs revealed in Theorem [31 

5. Convergence of ODEs' solution: a probabilistic approach 

Analogous to the steady-state probability distributions of the Markov 
chains underlying PEPA models, upon which performance measures such as 
throughput and utilisation can be derived, we expect the solution of the 
generated ODEs to have similar equilibrium conditions. In particular, if the 
solution has a limit as time goes to infinity we will be able to similarly obtain 
the performance from the steady state, i.e. the limit. Therefore, whether the 
solution of the derived ODEs converges becomes an important problem. 

We should point out that Kurtz's theorem cannot directly apply to the 
problem of whether or not the solution the derived ODEs converges. This 
is because Kurtz's theorem only deals with the approximation between the 
ODEs and Markov chains during any finite time, rather than considering the 
asymptotic behaviour of the ODEs as time goes to infinity. This section will 
present our investigation and results about this problem. 

5.1. Convergence under a particular condition 

We follow the assumptions in Theorem [3J Denote the expectation of 
X n (t) as M n (t), i.e. M n {t) = E[X n (t)]. For any t, the stochastic processes 
{X n (t)} n converge to the deterministic X(t) when n tends to infinity, as 
Theorem [3] shows. It is not surprising to see that {M n (t)} n , the expectations 
of {X n (t)} n , also converge to X(t) as n — > oo: 

Lemma 1. For any t, 

lim M n {t) =X(t). 

n— >oo 

Proof. Since X(t) is deterministic, then = X(t). By Theorem [31 

for all t, X n (t) converges to X(t) almost surely as n goes to infinity. Notice 
that X n (t) is bounded (see the discussion in the previous section), then by 
Lebesgue's dominant convergence theorem, we have 

lim E\\X n (t) - X(t)\\ = 0. 

n— too 
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Since a norm || • || can be considered as a convex function, by Jensen's in- 
equality (Theorem 2.2 in |3lj|). we have || (E[-]) \\ < E[\\ ■ \\]. Therefore, 

lim ||M n (*)-X(t)||= lim \\E[X n (t)]-E[X(t)]\\ 

n— yoo n— >oo 

< lim E\\X n (t)-X(t)\\ 

n—>oo 

= 0. 

□ 

Lemma [1] states that the ODEs' solution X(t) is just the limit function 
of the sequence of the expectation functions of the corresponding density 
dependent Markov chains. This provides some clues: the characteristics of 
the limit X{t) depend on the properties of {M n (t)} n . Therefore, we expect 
to be able to investigate X(t) by studying {M n (t)} n . 

Since M n (f) is the expectation of the Markov chain X n (t), M n (t) can 
be expressed by a formula in which the transient probability distribution is 
involved. That is, 



M n (t) = E[X n (t)} = ™ 



t 

x65„ 



where S n is the state space, 7r"(-) is the probability distribution of X n at 
time t. Let S n and tt?(-) be the state space and the probability distribution 
of X n (t) respect iveljo Then 



M n {t) = E[X n {t)} = E 



X n (t) 



n 



We have assumed the Markov chains underlying PEPA models to be irre- 
ducible and positive- recurrent. Then the transient probability distributions 
of these Markov chains will converge to the corresponding steady-state prob- 
ability distributions. We denote the steady-state probability distributions of 
X n (t) and X n {t) as 7r£,(-) and 7r^(-) respectively. Then, we have a lemma. 

Lemma 2. For any n, there exists a M n (oc), such that 

lim M n (t) = M„(oo). 

t— >oo 



2 We should point out that the probability distributions of X n (t) and X n (t) are the 
same, i.e. 7rf(x) = 7r"(x/n). 
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Proof. 

lim M n (t) 

t— >oo 

= lim -<w = E ; im -*?(*) = E -cw = ^n(oo). 

□ 

Clearly, we also have M n (oo) = X] x es n x ^oo( x )- 

Remark 2. Currently, we do not know whether the sequence {M n (oo)} n 
converges as n — > oo. -Efat£ since {M n (oo)} n is bounded which is due to the 
conservation law that PEPA models satisfy, there exists {n'} C {n} such that 
{M n /(oo)} converges to a limit, namely Moo(oo). That is 

lim M„/(oo) = Moo(oo). 

n'— too 

Thus, 

lim lim M n >{t) = lim M n ,(oo) = Moo(oo). (28) 

n'— >oo t— >oo n'— >oo 

At the moment, there are two questions: 

1. Whether lim^oo linv^^ M n i (t) exists? 

2. If lim^oo lim n /^.oo M n / (t) exists, whether 

lim lim M n >(t) = lim lim M n >{t)l 

t— >-oo n'— >oo n'— ¥oo t— >oo 

If the answer to the first question is yes, then the solution of the ODEs 
converges, since by Lemma [2} lim^oo X{t) = lim^oo linv^oo M n i(t). If the 
answer to the second question is yes, then the limit of X(t) is consistent with 
the stationary distributions of the Markov chains since 

lim X(t) = lim lim M n ,(t) = lim lim M n >(t) = Moo(oo). 

t— >oo t— >oo n'— >oc n'— >oo t— >oo 

In short, the positive answers to these two questions determine the conver- 
gence and consistency for the ODEs's solution, see Figure [U Fortunately, 
the two answers can be guaranteed by the condition ( 1291) in the following 
Proposition [6j 
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Figure 1: Convergence and consistency diagram for derived ODEs 



M n ,{t) 



Lemma [T] 



■» X(t)(=Moo(t)) 



Lemma [2] 



??? 



t— >0O 



M n ,(oo) 



Remark [2] 



Proposition 6. (A particular condition) If there exist A, B > 0, such 
that 



sup 



M n ,(t) - M n ,(oo) 



< Be 



-At 



(29) 



then lim 4 ^ 00 X(t) = M^oo) 



Proof. 



Xty-M^oo) 



lim [M n ,{t) - M n ,(oo)] 



< lim sup 



< lim sup 

n'— >oo 



M n ,(t)-M n ,(oc) 



sup 



M nf (t)-M n /(oo) 



< lim sup Be 

n'—toc 
-At 



Be' 



-At 



0, as t — > oo. 



So Hindoo X(£) = Moo(oo) 



□ 



Notice that MJt) = -7r, n (x) and M„(oo) = -7r" (x), so in order 



XG5" 



XGS" 



to estimate 



M n {t) — M n (oo) in (1291) . we need first to estimate the difference 



between 7r" and ir'^. 

Lemma 3. If there exists A > and i?i > 0, such that for any n' and all 
x G S n ' 

K'(x)-<(*)\<<(x)Bie- At , (30) 

< Be~ At holds. 



then there exists B > swc/i £/iai sup n / 



M n ,(t) -M n >(oo) 
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Proof. We know that X n (0) 



X n (0) 



n 



x for any n. By the conservation 



law, the population of each entity in any state is determined by the starting 
state. So for any n' and all xGS", ||x/n'|| < d J2pevM P ] 

< oo, where 

T> is the set of all local derivatives and C\ is a constant independent of n! . 
Let C = sup n / max xgS . n ' ||x/n'||, then C < oo. 



M n , (t) -M n , (<x>) 



X 

7l' 



xeS" 



xes™ 



< 



£ -(vrfW-TT^x)) 
z — ' n 

xes™' 

sup max ^ K'(x; 



xes™' 



xes™' 
xes™' 



< C > ' <(x)5 ie -^ 



Let 5 = CB 1 . Then sup n , M n /(t) - M n ,(oo 



< Be 



-At 



□ 



5.2. Investigation of the particular condition 

This section will present the study of the particular condition f )29|) . We 
will expose that the condition is related to well-known constants of Markov 
chains such as the spectral gap and the Log-Sobolev constant. The methods 
and results developed in the field of functional analysis of Markov chains are 
utilised to investigate the condition. 



5.2.1. Important estimation for Markov kernel 

We first give an estimation for the Markov kernel which is defined below. 
Let Q be the infinitesimal generator of a Markov chain X on a finite state 
S. Let 

1 + %, i = 3 



K l3 



where m = sup(— Qa). 
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K is a transition probability matrix, satisfying 

K(x,y)>0, £tf(x,y) = l. 
y 

K is called an Markov kernel (K is also called the uniformisation of the 
CTMC in some literature). A Markov chain on a finite state space S can be 
described through its kernel K. The continuous time semigroup associated 
with K is defined by 

H t = ex P (-t(I-K)). 

Let 7T be the unique stationary measure of the Markov chain. Then Ht(x, y) — > 
7r(y) as t tends to infinity. Following the convention in the literature we will 
also use (K, 7r) to represent a Markov chain. 
Notice 



Q = m{K-I), 



mt{K~I) 



-mt(I-K) 



Clearly, 

P t = e tQ = e mt ^- 1! = e - mt ^-«) = H mt , 
and thus H t = P±, where P t is called the semigroup associated with the 

m 

infinitesimal generator Q. An estimation of H t is given below. 

Lemma 4. (Corollary 2.2.6, fs^J) Let (K,ir) be a finite Markov chain, and 
7r(*) = min7r(x). Then 



sup 



# t (x,y) 



^(y) 



< e 



2-c 



for t = — log log 

a 



7T * 



c 

A' 



(31) 



where A > 0, a > are the spectral gap and the Log-Sobolev constant respec- 



tively, which are defined and interpreted in Appendix B 



It can be implied from fl3~Tl) that Vx, y G S 
|if*(x,y)-7r(y)| <7r(y) 
Pt , so 



2+^ log log -jL 

g a o o 7r(*) 



-A/ 



Since H, 



P±(x,y) - 7r(y) < 7r(y) 



2+^ lOg lOg -tr 

3 1 a O O 7r(*J 



-,\f 



and thus (replacing "t" by "mf on the both sides of (I33]T 



|P i (x,y)- 7 r(y)| <7r(y) 



2+^-loglog-f^ 

g 1 a O O 7r(*) 



-m\t 



(32) 

(33) 
(34) 



The investigation and utilisation of the formulae f[34|) will be presented in 
the next subsection. 
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5.2.2. Study of the particular condition 

For each n, let Q n be the infinitesimal generator of the density dependent 
Markov chain X n (t) underlying a given PEPA model and thus the transition 
probability matrix is P™ = e tc * n . For each X n (t), the initial state of the 
corresponding system is nx , so the initial probability distribution of X n (t) 
is Tr^(nxo) = 1 and <(x) = 0, Vx 6 S n , x ^ nx . So Vy G S n , 



xeS" 



<(x)Pf (x,y) = 7r >x )P t >x ,y) = P t >x ,y) 



The formula (1541) in the context of X n (t) is 



IA n (x,y)-^( y )|<^(y) 



2+^ log log n 1 , , 



-Tn n X n t 



(35) 



(36) 



where X n , X n , m n , 7r£>(*) are the respective parameters associated with X n (t). 

Let x = nx in ( 1561 . and notice P"(nx ,y) = 7r"(y) by ( 1551) . then we 
have 

2+^ log log n 1 , \ " 



W(y)-^(y)|<C(y) 



Vy e 5 n . (37) 



From the comparison of ( 1571) and (1501) . we know that if there are some 
conditions imposed on ^ log log T „^ and —m n X n t, then ( 150]) can be induced 
from (1571) . See the following Lemma. 

Lemma 5. // t/iere exists T>0,B2>0,A>0 such that 



sup 



-m n \ n T + — log log — 

Oirt. 



< Po 



and 
then 



inf{m n A n } > A, 



K(x)-7r»(x)|<C(x)B ie - 



Vx e 5" 



(38) 

(39) 
(40) 



where B\ = e A ^ T + B 2+ 2 ; an d ^ e "particular condition" 12~9\) holds. 
Proof. 



2+ An log l og 1 



-m„X n t 



-m^T+2+^loglog^y 



-m n X n (t-T) 



Pie 



-.1/ 



Then by (137)). (HO) holds. Thus by Lemma [5J d2HD holds. 



□ 
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Remark 3. When n tends to infinity, the discrete space is approximating a 
continuous space and thus the probability of any single point is tending to 7 
that is, 7r^(*) tends to 0. So log log oo as n goes to infinity. Notice 



that by LemmalT5\ in\Appendix B 



-<iog(iAC(*)). 

Therefore, in order to have 

-m n \ n T + — log log 1 < B 2 , 

«n 7T™(*) 

it is sufficient to let 

Tm n A n >0([log(l/^(*))] 2 ). (41) 
Moreover, can imply both [33^1 and / fffPj) . 

According to the above analysis, our problem is simplified to checking that 
whether (14T|) is satisfied by the density dependent Markov chains {X n (t)}. 
By Remark [H] in Appendix B X n is the smallest non-zero eigenvalue of 

j _ K n +K*' n ^ ^gj-g 

O n 

K n = — h I 

and K*' n is adjoint to K n . A matrix Q* ,n is said to be adjoint to the generator 
Q n , if g n (x,y)7r^(x) equals Q*' n (y, x)vr™ (y). Clearly, Q*> n = m n (K*> n - /), 
or equivalently, 

/n*,n 

K*' n = ^— + /. 
m n 

So 

m»(/ g )= g • (42) 

Denote the smallest non-zero eig envalue of - °" + 2 Q *' m by a n . Then by (H2jl . 

m n A n = <r n . (43) 
Now we state our main result in this section. 
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Theorem 4. Let {X n (t)} be the density dependent Markov chain derived 
from a given PEPA model. For each n G N, let S n and ir^ be the state 
space and steady-state probability distribution of X n (t) respectively. Q n is the 
infinitesimal generator of X n (t) and a n is the smallest non-zero eigenvalue 
of — ® + Q ' , where Q* ,n is adjoint to Q n in terms of ir^. If 



min 7r" (x) > 



exp(0(y/d~n)) 



(44) 



for sufficiently large n, then X(t) has a finite limit as time tends to infinity, 
where X(t) is the solution of the corresponding derived ODEs from the same 
PEPA model. 

Proof. By the given condition of vr™ (*) > exp(0 ^ )} , 

1 



log 



Thus 



< log[exp(0(V^))] = 0{a]! 2 ). 



< OK). 





1 


^log 









Choose a large T such that 





1 


flog 






.^(*). 



By Remark [3] and Lemma 0, the particular condition holds. Therefore 



lim X(t) 

t—too 



M^(oo) 



□ 



According to the above theorem, our problem is simplified to check- 
ing that whether (I44p is satisfied by the density dependent Markov chains 
{X n (t)}. In (jUJ) both the spectral gap A n (= — ) and 7r£>(*) are unknown. 
In fact, due to the state space explosion problem, 7r^ cannot be easily solved 
from n^Q 12 = or equivalently ir^ >0 K n = ir^. Moreover, the estimation of 
the spectral gap for a given Markov chain in current literature (e.g. 32]) is 



heavily based on the known stationary distribution. Thus, the current results 
cannot provide a practical check for 
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The convergence and consistency are supported by many numerical ex- 
periments, so we believe that ( I44p is unnecessary. In other words, we believe 
that ()44p always holds in the context of PEPA, although at this moment we 
cannot prove it. Before concluding this section, we leave an open problem: 

Conjecture 1. The formula [J4\ ) in Theorem^ holds for any PEPA model. 



6. Fluid analysis: an analytic approach (I) 

The previous sections have demonstrated the fluid approximation and rel- 
evant analysis for PEPA. Some fundamental results about the derived ODEs 
such as the boundedness, nonnegativeness and convergence of the solutions, 
have been established through a probabilistic approach. In this section we 
will discuss the boundedness and nonnegativeness again, and prove them by 
a purely analytical argument. The convergence presented in the previous 
section is proved under a particular condition that cannot currently be eas- 
ily checked. This section will present alternative approaches to deal with 
the convergence problem. In particular, for an interesting; model with two 



synchronisations, its structural invariance as revealed in [22j, will be shown 
to play an important role in the proof of the convergence. Moreover, for a 
class of PEPA models which have two component types and one synchroni- 
sation, an analytical proof of the convergence under some mild conditions on 
the populations will be presented. These discussions and investigations will 
provide a new insight into the fluid approximation of PEPA. 

6.1. Analytical proof of boundedness and nonnegativeness 

Recall that the set of derived ODEs from a general PEPA model is 

^ = - E fM+ E /(*•<>• (45) 
0|t/epre(0} {i\uepost(i)} 

As mentioned, in this formula the term ^{;|^ e pre(z)} /( x ' represents the 
exit rates in the local derivative U . An important fact to note is: the exit 
rates in a local derivative Cj . in state x are bounded by all the apparent rates 
in this local derivative. In fact, according to Proposition [1] in Section \2?2[ if 
Cj. G pre(Z) where I is a labelled activity, then the transition rate function 
/(x, /) is bounded by r/(x, C^), the apparent rates of I in Cj. in state x. That 
is, 

/(x,Z)<r f (x,(7 i .)=x[Ci i h(^). (46) 
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We should point out that ( l4"6"j) is based on the discrete state space underlying 
the given model. According to our semantics of mapping PEPA model to 
ODEs, the fluid approximation-version of (I46p also holds, i.e. /(x(t),Z) < 
x(CV, t)ri(Ci ). Hereafter the notation x[-] indicates a discrete state x, while 
x(-) or x(-,t) reflects a continuous state x at time t. Therefore, we have the 
following 

Proposition 7. For any local derivative Ci j; 

/(x(t),Z)<x(a.,t) £ n(C y, (47) 
{i\Ci.epre(i)} {i\Ci.epre(i)} 

where 77 (Cj.) is the apparent rate of I in for a single instance of C ir 
defined in Definition^ 

Proposition [7] and Propositions [3] can guarantee the boundedness and non- 
negativeness of the solutions. In the following, we will present an analytical 
proof of these properties, based on the two propositions. Suppose the initial 
values x [C^ , 0) are given, and we denote ^ ■ x [Ci. , 0) by Nc 4 ■ We have a 
theorem: 

Theorem 5. If x (Cj. , t) satisfies ( f75| ) with nonnegative initial values, then 

< x (d, , t) < N Ci , for any t > 0. (48) 

Moreover, if the initial values are positive, then the solutions are always 
positive, i.e., 

< x (d, , t) < N Ci , for any t > 0. (49) 

Proof. By Proposition [3j ^2,- x (pi- , tj = Nd for all t. All that is left to do 
is to prove that x(Cj ., t) is positive or nonnegative. The proof is divided into 
two cases. 

Case 1: Suppose all the initial values are positive, i.e. min^ |x(Cj., 0)} > 
0. We will show that min^. {x(C^,t)} > for all t > 0. Otherwise, if there 
exists a t > such that min^. |x(CV,.,t)} < 0, then there exists a point t' > 
such that min^. |x(Ci j ,t / )} = 0. Let t* be the first such point, i.e. 

t* = inf i t > | min {x(C ij , t)} = i , 
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then < t* < oo. Without loss of generality, we assume x(Ci 1 ,t) reaches 
zero at t*, i.e., 

x(c ll ,r) = o, x(a„O>0 (^ivjvi) 

and 

x (a J ,t)>o, tG[o,r), vi,j. 

Thus, for t G [0,t*], by Proposition U\ 

= - £ /< x .')+ £ /(*.') 

{«|Ci ie pre(0} oic^epostcz)} 

> - £ /( x ' z ) 

{i|cr Xl epre(0} 

> -x(C ll5 t) £ rKCiJ. 

{«|c ll6 pre(/)} 

Set i2 = E{i|«7i ie pre(0} r «(di), then 



> -ifc^d^i). (50) 



By Lemma [T4l in Appendix A (l50l) implies 

x(C ll( t*) >x(C ll ,0)e^ i * >0. 

This is a contradiction to x(Ci 1 ,t*) = 0. Therefore < x (C^. and thus 
by Proposition [31 

0<x(C tj ,f) <N Co , Vt. 

Case 2: Suppose min^ {x(Cj j , 0)} = 0. Let w^i,-, 0) = x(Cj j ., 0) + 6 where 
5 > 0. Let us(ij,t) be the solution of ( )45|) . given the initial value ug(ij,0). 
By the proof of Case 1, us(ij,t) > (V£ > 0). Noticing min(-) is a Lipschitz 
function, by the Fundamental Inequality in Appendix A, we have 

\u s (ij,t)-yi(C ij7 t)\<8e K \ (51) 

where K is a constant. Thus, for any given t > 0, 

x(C ij ,t) > u*(ij,t) - <?e** > -<5e^. (52) 

Let 5 4- in (152|) . then we have x(C» ,, t) > 0. The proof is completed. □ 
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6.2. A case study on convergence with two synchronisations 

If a model has synchronisations, then the derived ODEs are nonlinear. 
The nonlinearity results in the complexity of the dynamic behaviour of fluid 
approximations. However, for some special models, we can still determine the 
convergence of the solutions. What follows is a case study for an interesting 
PEPA model, in which the structural property of invariance will be shown 
to play an important role in the proof of the convergence. 

6.2.1. Theoretical study of convergence 

The model considered here is given below, 



Model 2. 



x 2 

Yi 

Y s 
Y 3 
Y, 



{X^MAWX^Mg]) 



def 
def 
def 
def 
def 
def 

action! ,action2 



(action 1 , a{).X 2 
(action2, a 2 ).X 1 

(action 1 , ai). Y 3 + (jobl , cj). Y 2 

(job2,c s ).Yt 

(job3,c 3 ).Y 4 

(action2, a 2 ).Y 2 + (job4, c 4 ).Y 3 
(F,[^]||r 2 [iV 2 ]||r s [iV s ]||^[i^]) 





Figure 2: Transition systems of the components of Model 

The operations of X and Y are illustrated in Figure [2j According to the 
mapping semantics, the derived ODEs from this model are 
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where 



"dt" 

dx 2 

~~dT 

dt 
dy 2 
dt 
dy 3 
dt 

dt 
1,2, 



— ai min{xi, + a 2 min{a;2, t/ 4 } 
ai min{xi, y{\ — a 2 min{x2, 7/4} 
-ai min{a;i, J/i} - cij/i + c 2 y 2 
-c 2 y 2 + a 2 min{x 2 , t/ 4 } + c x y\ 
-c 3 y-3 + c 4 t/ 4 + ai min{a;i, 2/1} 
-a 2 min{x 2 , 7/4} - c 4 ?/4 + c 3 y 3 



(53) 



J 



1,2, •■■ ,4) denote the populations of X and 
Y in the local derivatives Yj respectively. Now we state an interesting 
assertion for the specific PEPA model: the difference between the number 
of Y in their local derivatives Y3 and Y4, and the number of X in the local 
derivative X 2 , i.e. y 3 + y& — x 2 , is a constant in any state. This fact can be 
explained as follows. Notice that there is only one way to increase 2/3 + 7/4, i.e. 
enabling the activity actionl. As long as actionl is activated, then there is a 
copy of Y entering the Y 3 from Y\. Meanwhile, since actionl is shared by X, a 
corresponding copy of X will go to X 2 from X\. In other words, 7/3+7/4 and x 2 
increase equally and simultaneously. On the other hand, there is also only one 
way to decrease 7/3 + 7/4 and x 2 , i.e. enabling the cooperated activity action2. 
This also allows 7/3 + 7/4 and x 2 to decrease both equally and simultaneously. 
So, the difference 7/3+7/4— x 2 will remain constant in any state and thus at any 
time. The assertion indicates that each state and therefore the whole state 
space of underlying CTMC may have some interesting structure properties, 
such as invariants. The techniques and applications of structural analysis of 



PEPA models have been developed in [22 



Throughout this section, let M and be the total populations of the 
X and Y respectively, i.e. M = M 1 + M 2 and N = Aq + N 2 + N 3 + JV 4 . 
Notice 7/1 + 7/2 + 7/3 + 7/4 = iV and x\ + x 2 = M by the conservation law, so 
7/1 + 7/2 — x\ = N — M — (7/3 + 7/4 — x 2 ) is another invariant because 7/3 + 7/4 — x 2 
is a constant. The fluid- approximation version of these two invariants also 
holds, which is illustrated by the following Lemma [6j 
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Lemma 6. For any t >0, 



yi(t) + y 2 (t) - x^t) = yi (o) + y 2 (o) - Xi(o), 
y 3 (t) + y 4 (f) - x 2 (t) = 2/ 3 (0) + 2/4(0) - x 2 (0). 

Proof. According to (|53|) . for any t > 0, 

d(yi(t) +V2(t) -xi(t)) _ dyi dy 2 dxi _ 
dt dt dt dt 

So 2/1 (t) H- 2/2 (*) — xi(t) = 2/i(0) +2/2(0) —£1(0), Vt > 0. By a similar argument, 
we also have 2/ 3 (t) + y 4 (t) - x 2 (t) = 2/3(0) + 2/4(0) - x 2 (0), Vt > 0. □ 

In the following we show how to use this kind of invariance to prove the 
convergence of the solution of fl53l) as time goes to infinity Before presenting 
the results and the proof, we first rewrite (|53|) as follows: 



/ dgi \ 
1 dt \ 

dx2 
dt 
ayi 
dt 
dy2 

ddt 

am 

dt 

\ ^y± I 

\ dt ' 



l {xi<y 1 ,x 2 <yi}Ql 



( X l \ 

X 2 

2/i 

2/2 
2/3 

V 2/4 / 
fx, > 

•''2 

2/1 

2/2 
2/3 

V 2/4 y 



! {^i<yi^2>j/4> ( 32 



+ ^{xi>2/i,a;2>2/4}Q4 



/ Xi \ 

X 2 
tJl 
2/2 
2/3 

V 2/4 / 

/ Xi \ 

x 2 
2/i 
2/2 

2/3 

V 2/4 / 



(54) 



where the matrices Qi (i = 1, 2, 3, 4) are given as below: 



Qi 



—d\ 


a-2 


\ 


a. 


—a 2 




—a\ 




-Cl c 2 




a 2 


Cl -c 2 






— C3 C4 




-a 2 


c 3 -c 4 / 
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g 2 = 



/ 


— d\ 









a 2 


\ 




ai 









— a 2 






-ai 





-Cl c 2 















Cl -c 2 




a 2 






ai 







-c 3 


c 4 




V 









c 3 


-(c 4 + a 2 ) 


/ 



/ o 

_0_ 







«2 

-a 2 
a 2 
-a 2 



— ai 
ai 
-ai — Ci 

Cl 



c 2 
-c 2 



-c 3 

c 3 



c 4 
-c 4 / 



Ql 



( 





— ai 








ai 








— ai — ci 








Cl 








ai 










c 2 
-c 2 



-c 3 

c 3 



a 2 
— a 2 

a 2 
c 4 

-a 2 - c 4 / 



As (1541) illustrates, the derived ODEs are piecewise linear and they may 
be dominated by Qi (i = 1, 2, 3, 4) alternately. If the system is always domi- 
nated by only one specific matrix after a time, then the ODEs become linear 
after this time. For linear ODEs, as long as the eigenvalues of their coefficient 
matrices are either zeros or have negative real parts, then bounded solutions 
will converge as time tends to infinity, see Corollary [3] in Appendix D For- 
tunately, here the eigenvalues of the matrices Qi {i = 1,2,3,4) in Model [2] 
satisfy this property, the proof of which is shown in Appendix C| In addi- 
tion, the solution of the derived ODEs from any PEPA model is bounded, 
as Theorem \5\ illustrated. Therefore, if we can guarantee that after a time 
the ODEs f!54p become linear, which means that one of the four matrices 
Qi (i = 1,2,3,4) will be the coefficient matrix of the linear ODEs, then by 
Corollary [3] the solution will converge. So the convergence problem is reduced 
to determining whether the linearity can be finally guaranteed. 

It is easy to see that the comparisons between x\ and yi, x 2 and y 4 de- 
termine the linearity. For instance, if after a time T, we always have X\ > y\ 
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and %2 > 1/4, then the matrix Q 4 will dominate the system. Fortunately, 
the invariance in the model, as Lemma E] reveals, can determine the com- 
parisons in some circumstances. This is because this invariance reflects the 
relationship between different component types that are connected through 
synchronisations. This leads to several conclusions as follows. 

Proposition 8. Ifyi{0) + 2/2(0) < £i(0) and y 3 (0) + 2/4(0) < x 2 (0), then the 
solution of converges. 

Proof. By Lemma [61 yi(t) + 2/2 (*) < xi(t) and 2/3 (t) + y±(t) < x 2 (t) for all 
time t. Since both 2/ 2 (t) and y&{f) are nonnegative by Theorem |5l we have 
yi(t) < xi(t) and 2/4 (t) < X2(t) for any t. Thus, fl54"j) becomes 

/dxi dx 2 dj/i dy 2 dy 3 dyA T T 
(.IT'lT' dP dT' dT' dTj =^ (^,^,^,^,2/3,2/4) . (55) 

Notice that ( 155|) is linear, and all eigenvalues of Q4 other than zeros have neg- 
ative real parts, then according to Corollary [3j the solution of (|55|) converges 
as time goes to infinity. □ 

Proposition 9. Suppose y\(0) +2/2(0) > xi(0) and 2/3(0) +2/4(0) < x 2 (0). 
// either (I). N > ( 2 + !H±?L) M, or (II). N > 2 ^ + c ^ + a * Mj where 

V C 2 / c 2 

N > M > are £/ie populations of Y and X respectively, then the solution 
of |5^P converges. 

Proof. Suppose (I) holds. According to the conservation law, Y^i=i Vi{t) = N. 
By the boundedness of the solution, we have x 2 (t) < M. Then by Lemma El 
ya(*) + 2/4(t) < x 2 (t) < M. Therefore, 

y 2 (t) = N- (y 3 (t) + y 4 (t) ) — yi{t) > N — M — Vl (t). (56) 
Since min{xi,2/i} < 2/1, so — a min{xi, 2/1} > — ai2/i- Thus 

— = -ai mmjxi, 2/1} - C12/1 + c 2 y 2 

> -am - cij/i + c 2 2/ 2 (57) 

> -(ai + ci)2/i + c 2 (iV - M - 2/0 
= -(ax + ci + c 2 )2/i + c 2 (iV - M). 
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That is 

^ > -(ai + ci + c 2 ) yi + c 2 (N - M). 
Applying Lemma [TH in Appendix A| to this formula, we have 

c 2 (N-M) V_ (oi+Cl+C2)t , c 2 (N-M) 



Vl (t) > Vl (0) - ^\ ' e-^*** + \ (58) 

\ ai + ci + c 2 / ai + Ci + c 2 

Since the first term of the right side of (1581) converges to zero as time goes 

to infinity, i.e. lim ( yAO) - ° 2 ^ N ~ M > ) e -(«»i+ci+Q9)t = q, and the second 
V ai + ci + c 2 / 

term — ^ > M which results from the condition N > \ 2 -\ — ^— — - ] M, 

ai + ci + c 2 V c 2 / 

then there exists T > such that for any t > T, y\(t) > M > Xi(t). Then 
after time T, 054|) becomes linear, and is dominated by Q 2 . Because all 
eigenvalues of Q 2 are either zeros or have negative real parts, the solution 
converges as time goes to infinity. 

Now we assume (II) holds. Similarly, since min{x 2 ,?/4} < x 2 < M, and 
yi < N — y 2 which is due to y 1 + y 2 < N, we have 

-7T = - c 2Z/2 + a 2 mm{x 2 , y 4 } + cij/i 
d£ 

< -c 2 y 2 + a 2 M + cij/i (59) 

< -c 2 y 2 + a 2 M + Cl (A-y 2 ) 
= -(ci + c 2 )y 2 + a 2 M + ciJV. 



By Lemma in Appendix A 



/ a 2 M + ciA\ _ (ci+C2) , a 2 M + ciA^ 
2/2 < 1/2(0) ■ e (C1+C2)i + ■ . (60) 

V ci + c 2 y ci + c 2 

Therefore, since e~^ cl+C2 ^ in above formula converges to zero as time tends 
to infinity, then for any e > 0, there exists T > such that for any time 
t > T 

a 2 M + Cl N 

V2 < ■ + e. 61 

ci + c 2 

Notice that the condition N > - 1 2 ^ — -M implies 

c 2 

c 2 N - a 2 M - (ci + c 2 )M M 
Cl + c 2 
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hi 11 11 c iN - a 2 M - (ci + c 9 )M , , 

and let e be small enough that e > M. lhen by 

ci + c 2 

§B, Vi> (N - M) -y 2 . Therefore, 

yi >(N- M) - y 2 

>(iV-M)- a2M + CliV - £ 

Cl + ° 2 (62) 
_ c 2 N - a 2 M - (ci + c 2 )M 1 ; 

Cl + c 2 

> M > xx. 



So 2/1 (t) > X\(t), 7/4 (t) < X2(t), for any t > T, then by a similar argument 
the solution of (154D converges. □ 



Both condition (I) and (II) in Proposition require to be larger enough 
than M, to guarantee that y\ is larger than x\. Since our PEPA model is 
symmetric, Proposition |9] has a corresponding symmetric version. 

Proposition 10. Suppose ?/i(0) + 2/2(0) < xi(0) and 2/3(0) + 2/4(0) > a; 2 (0). 

// either (I). N > ( 2 + M ; or (71). iV > 2( ° 3 + ° a) + Ql M ; ^ere 

V c 4 / Ci 

A 7 " > M > are the populations of Y and X respectively, then the solution 
of ( [5^] ) converges. 

The proof of Proposition [10] is omitted here. We should point out that in 
our model the shared activity actionl (respectively, action2) has the same 
local rate a\ (respectively, a 2 ). We have taken the advantage of this in the 
above proofs. If the local rates of shared activities are not set to be the 
same, analogous conclusions can still hold but the discussion will be more 
complicated. However, the structural property invariance can still play an 
important role. 

The above three propositions have illustrated the convergence for all sit- 
uations in terms of the starting state, except for the case of 2/1 (0) + 2/2(0) > 
^1(0), 2/3(0) + 2/4(0) > x 2 (0). See a summary in Table [3j If 2/i(0) + 2/2(0) > 
x i(0), 2/3(0) + 2/4(0) > £2(0), then the dynamic behaviour of the system is 
rather complex. A numerical study for this case will be presented in the next 
subsection. 
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Table 3: A summary for the convergence of Model [5] 



Starting state condition 


Additional condition 


Conclusion 


yi(0)+y 2 (0) < Xl (0), 
|/3(0)+ 2/ 4 (0) <x 2 (0) 


No 


Proposition |8] 


2/i(0) + 2/2(0) >x 1 (0), 
2/3(0) + 2/4(0) <x 2 (0) 


iV > fciM 


Proposition |9] 


2/i(0) +2/2(0) <2i(0), 

2/3(0) + 2/4(0) >x 2 (0) 


N > k 2 M 


Proposition [10] 


2/i(0) + 2/ 2 (0) >zi(0), 
2/3(0) + 2/4(0) >x 2 (0) 


None identified 


Explored numerically 




0123456789 10 

time 



Figure 3: Numerical study for Model [5] rates (1,1,1,1,1,1); equilibrium point 
(1, 1, 2, 3, 3, 2). (Note: the curves of x\ and X2, the curves of y\ and 2/4, as well as those of 
1/2 and 2/3 respectively, completely overlap.) 
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time 



Figure 4: Numerical study for Model H) rates (1,1,1,10,1,10); equilibrium point 
(0.4616, 1.5384, 4.0140, 0.4476, 5.0769, 0.4615) 




time 



Figure 5: Numerical study for Model [2j rates (1,1,10,1,1,1); equilibrium point 
(1.5384, 0.4616, 0.4615, 5.0769, 2.4616, 2.0000) 
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5 
time 



Figure 6: Numerical study for Model [2] rates (20,20,1,1,1,1); equilibrium point 
(1,1,0.2273,4.7727,4.7727,0.2273). (Note: the curves of xi and x 2 , the curves of yt 
and ?/4, as well as those of yi and j/3 respectively, completely overlap.) 



6.2.2. Experimental study of convergence 

This subsection presents a numerical study at different action rate condi- 
tions. The starting state in this subsection is always assumed as (1, 1, 5, 0, 0, 5), 
which satisfies the condition of 7/1 (0) +1/2(0) > £i(0) and 1/3(0) +7/4(0) > a; 2 (0). 

If all the action rates in the model are set to one, i.e. (ai, a 2 , c\, c 2 , c 3 , c 4 ) = 
(1, 1, 1, 1, 1, 1), then the equilibrium point of the ODEs is (x{, x* 2 , y\, y%, y%,yl) = 
(1, 1, 2, 3, 3, 2), as the numerical solution of the ODEs illustrates. In this case, 
the matrix Qi finally dominates the system. See Figure |3j Notice that in 
this figure, the curves of x\ and x 2 completely overlap, as well as the curves 
of 2/1 and 7/4, and those of y 2 and y 3 . 

In other situations, for example, if (a l5 a 2 , ci, c 2 , c 3 , c 4 ) = (1, 1, 1, 10, 1, 10) 
then the equilibrium point is (0.4616,1.5384,4.0140,0.4476,5.0769,0.4615), 
and the matrix Q 2 eventually dominates the system. See Figure HJ Moreover, 
the matrices Q3 and Q 4 can also finally dominate the system as long as the 
action rates are appropriately specified. See Figure [5] and Figure |6j We 
should point out that similarly to Figure [31 in Figure M the curves of X\ 
and x 2 , the curves of y\ and y 4 , as well as those of y 2 and 773 respectively, 
completely overlap. 
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In short, the system dynamics is rather complex in the situation of yi(0) + 
2/2(0) > £i(0) and 1/3(0) + 2/4(0) > x 2 (0). A summary of these numerical 
studies is organised in Table HJ 



Table 4: Complex dynamical behaviour of Model [5J starting state (1,1,5,0,0,5) 



Rates: 
(ai,a2,ci,C2,c 3 ,c 4 ) 


Equilibrium points: 

/ S|S ^ 5f! + + + \ 

^1 , x 2 > 2/l > 2/2 > J/3 > 2/4 J 


Dominator 


Figure 


(1,1,1,1,1,1) 


(1,1,2,3,3,2) 


Qi 


Figure [3] 


(1,1,1,10,1,10) 




Q2 


Figure H] 


(1,1,10,1,1,1) 




Q3 


Figure [5] 


(20,20,1,1,1,1) 


(1, 1, 0.2273, 4.7727, 4.7727, 0.2273) 


Qa 


Figure M 



a vi = (0.4616,1.5384,4.0140,0.4476,5.0769,0.4615) 
% = (1.5384,0.4616,0.4615,5.0769,2.4616,2.0000) 



7. Fluid analysis: an analytic approach (II) 

7.1. Convergence for two component types and one synchronisation (I): a 
special case 

The problem of convergence for more general models without strict con- 
ditions, is rather complex and has not been completely solved. But for a par- 
ticular class of PEPA model — a model composed of two types of component 
with one synchronisation between them, we can determine the convergence 
of the solutions of the derived ODEs. 

As discussed in the previous section, the ODEs derived from PEPA are 
piecewise linear and may be dominated by different coefficient matrices al- 
ternately. For any PEPA model which has two component types and one 
synchronisation, the two corresponding coefficient matrices can be proved 
to have a good property: their eigenvalues are either zeros or have negative 
real parts. The remaining issue for convergence is to ascertain that these 
two matrices will not always alternately dominate the system. In fact, we 
will prove that under some mild conditions, there exists a time after which 
there is only one coefficient matrix dominating the system. This means the 
ODEs become linear after that time. Since the coefficient matrix of the lin- 
ear ODEs satisfies the good eigenvalue property, then by Corollary El the 
bounded solution will converge as time goes to infinity. 
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We first utilise an example in this section to show our approach to dealing 
with the convergence problem for this class of PEPA models. The proof for 
a general case in this class is presented in the next subsection. 

7.1.1. A previous model and the fluid approximation 

Let us look at the following PEPA model, which is Model [T] presented 
previously: 

Useri =(taski,a).User 2 

User 2 = {tash 2 ,b).User 1 
Provider \ = (task\, a). Provider 2 
Provider 2 = (reset,d).Provideri 
[User^M]) D*a (Provider^N}) . 

The derived ODEs follows: 

dxi 



dt 
dx 2 

~dT 

dt 
dy 2 
dt 



—a min{xi, y±} + bx 2 
a min{xi, y{\ — bx 2 
-amm{x ll y 1 } + dy 2 
amin{xi,?/i} - dy 2 



(63) 



where Xi and yi represent the populations of Useri and Provideri respec- 
tively, i = 1,2. Clearly, (!63j) is equivalent to 



where 



/ d£i \ 

/ dt 1 

dx2 
dt 

am 
dt 

\ zuz I 

\ dt / 

/ 



I{x 1 <y 1 }Ql 



( X l \ 

yi 

\V2 ) 



Qi 



—a b 





\ 


a —b 








-a 





d 


a 





-d) 



( b 
-b 




\ 



( Xi \ 

x 2 

yi 

\ yi ) 



-a \ 

a 



(64) 



—a d 
a —d ) 



(65) 
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Our interest is to see if the solution of (I64p will converge as time goes to 
infinity. As we mentioned, this convergence problem can be divided into two 
subproblems, i.e. whether the nonlinear equations can finally become linear 
and whether the eigenvalues of the coefficient matrix are either zeros or have 
negative real parts. If answers to these two subproblems are both positive, 
then the convergence will hold. 

The second subproblem can be easily dealt with. By calculations, the 
matrix Qi has eigenvalues (two folds), —d, and —(a + b). Similarly, Q 2 
has eigenvalues (two folds), —6, — (a + d). Therefore, the eigenvalues of Qi 
and Q2 other than zeros are negative. Moreover, for a general PEPA model 
which has two component types and one synchronisation, Theorem [7] in the 
next section shows that the corresponding coefficient matrices always have 
this property. 

The remaining work to determine the convergence of the ODE solution, 
is to solve the first subproblem, i.e. to ascertain that after a time it is always 
the case that X\ > y\ or X\ < y±. In this model, there is no invariance 
relating the two different component types, so we cannot rely on invariants 
to investigate this subproblem. However, we have a new way to cope with 
this problem. 

7.1.2. Proof outline of convergence 

Notice that y\{t) < N by the boundedness of solutions. If we can prove 
that after time T, Xi(t) > cM, where c > is independent of M, we will get, 
provided cM > N, 



Therefore, the ODEs fl63|) will become linear after time T. In the following, 
we identify that Xi(t) > cM. 



then < a(t) < 1 by the nonnegativity of Xi(t) and yi(t). The ODEs 
associated with component type X can be rewritten as 



xi{t) >cM>N> yx(t), t > T. 



Let 




xi (t) ^0, 
x 1 (t) = 0, 




(66) 
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Let 



w - 1 It -6 ) • ^ - ( It 



Then (1661) can be written as 
The solution of (1671) is 



A(t)X(t). (67) 



X(t) =e^^) ds X(0). (68) 



1 " 



Let 5(t) = - / A(s)ds, then 



o 



X(t) = eti A{s)ds X(0) = e tB{t) X(0), (69) 

and 



, . If -af Q a(s)ds bt \ _ ( -a(3{t) b \ 

1 j ~ t V afia(s)ds -bt J~\ a(3(t) -b )> [m) 

, , f* a(s)ds 

where /3(t) = Jo - — . Obviously, < (5{t) < 1 because < a{s) < 1. 
Notice that the matrix -B(t) can be diagonalised as 

B(f) = 



-a(3(t) b 
a/3(t) -b 

b i 



\ 1 1 



a %t fe _i I ( o -fa^m + M J I am — J ^ 



where 



6 1 \ II 1 



=( a % h , IT 1 ® = 6_ , (72) 

\ aj9(i)+& 1 / \ a/3{t)+b af3{t)+b J 

SO 

= e^WX(O) = 17(*) ( I e - t{a ° m+b) ) U~\t)X(0). (73) 
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In the formula (17TI) . both and — (af3(t) + b) are -B(t)'s eigenvalues, 
while the columns of U(t) are the corresponding eigenvectors. In particular, 

rp 

i.e. the first column of U(t), is the eigenvector 



a/3{t) + b' a/3{t) + b 
corresponding to the eigenvalue zero. 
We define a function X(t) by 



X(t) = U(t)(l ° ) U-\t)X(0). (74) 

By simple calculation, 

~ , , , 4M , M , T / bM aB(t)M \ T 

Clearly, is the normalised eigenvector corresponding the zero eigenvalue 
(at the time t). Here the normalisaton is in terms of the total population 
M of the component type X. Moreover, X(t) embodies some equilibrium 
meaning. In fact, we have a conclusion: 

\im\\X(t)-X(t)\\=0. (76) 

t— >oo 

Since the explicit expression of U(t) is available, the proof of (j76l) is easy 
and thus omitted. This conclusion is also included in Proposition [121 the 
proof of which does not rely on the explicit information of U(t) and will be 
introduced later. 

Now we discuss the benefit brought by this formula. By (ITS"]) the first 

entry of X(t) approximates the first entry of X(t), i.e. X\(t) approximates 
bM 

£i(t) = m 7- Thus, for any e > 0, there exists T > such that for any 

ap(t) + b 

t > T, 

Xl{t)>Xlit) - e = ^{tHb- e - 

bM bM W1 , . bM 

Since < — — < M because < (3(t) < 1, so xAt) > - e. 

a + b ap(t) + o a + b 

Therefore, if > N, then by the boundedness of yi(t), i.e. yi(t) < N, 

a + b 

we have 

Xl * >—r-e>N>yi(t), 
a + b 
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as long as e is small enough. This means that Q2 will dominate the system 
after time T. So we have 



Proposition 11. If 



bM 
a + b 



> N, then the solution of the ODEs converges 



as time tends to infinity. 

Since the model is symmetric, this proposition has a symmetric version: 



As we discussed, there are two key steps in the proof of Proposition [TT1 
The first step is to establish the approximation of xx(t) to Xi(t), i.e. X\{i) ~ 
Xi(t). The second one is to give an estimation Xi(t) > cM. According to 
these two conclusions, we have Xi(t) > c'M where d < c, and therefore can 
conclude that X\(t) > dM > N > y\(t) provided the condition c'M > N. 
This is the main philosophy behind our proof for x\(t) > yi(t). 

For the sake of generality, the proofs of these two conclusions should not 
rely on the explicit expressions of the eigenvalues and eigenvectors of B(t). 
This is because for general PEPA models with two component types and 
one synchronisation, these explicit expressions are not always available. The 
following subsection will present our discussions about these steps, and the 
proofs for the conclusions which do not rely on these explicit expressions. 

7.1.3. Proof not relying on explicit expressions 

This subsection will divide into two parts. In the first part, we will give 
a lower bound for the eigenvalues of the coefficient matrix B(t), based on 
which a proof of the approximation of X(t) to X(t) is given. The second 
part will establish the estimation X\(t) > cM. All proofs in this subsection 
do not require knowledge of the explicit expressions of the eigenvalues and 
eigenvectors of B(t). 

For convenience, in this subsection we define 



where / is a matrix-valued function defined on R. Then the matrix 



if 



dN 
a + d 



> M, then the solution of the ODEs also converges. 




(77) 
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can be written as B(t) = f((3(t)). The diagonalisation of f{(3) is 



where A(/3) is /(/3)'s nonzero eigenvalue, and g((3) is a matrix whose columns 
are the eigenvectors of f{(3). Here g~ l {[3) is the inverse of the matrix g(f3). 
Notice that A(/3) is real, because if A(/3) is complex then its conjugation must 
be an eigenvalue, which is contradicted by the fact that f{(3) only has two 
eigenvalues, and A(/3). The following discussions in this subsection do not 
rely on the explicit expressions of A(/3), g{(3) and g _1 (/3), although it is easy 
to see that A(/3) = —a(3 + b and 

m = ( 3? \ ), g-\P) = ( ]§_ ) ■ 

\ a/3+b 1 J \ afi+b a/3+b J 

1. X(t) approximates X(t) 



In the following, we will give a lower bound for the nonzero eigenvalue of 
f(P), i.e. A(/3), and based on this prove the approximation of A(t) to X(t) 
as time tends to infinity. 

If P > 0, then the transpose of /(/?), i.e. f((3) T , is an infinitesimal 
generator, and thus the nonzero eigenvalue A(/3) has negative real part, i.e. 
!R(A(/3)) < 0. If (3 = 0, then f((3) is independent of (3 and becomes a 
nonnegative matrix, i.e. each entry of it is nonnegative. Based on the Perron- 
Frobenious theorem which is presented in the next subsection, we can still 
have 9ft(A(0)) < 0. Therefore, for any j3, /(/3)'s eigenvalue other than zero 
has negative real part. This conclusion is stated in the following lemma. 

Lemma 7. For any (3 £ [0, 1], 9ft(A(/3)) < 0, where A(/3) is a nonzero eigen- 
value of f((3). 



The proof of Lemma [7] is presented in Appendix D.3 Lemma [7] can 
further lead to the following 

Lemma 8. Let 

Ai = inf {— 9ft(A(/3)) | X((3) is f((3) 's non-zero eigenvalue} , (78) 
/3e[o,i] 

then Ai > 0. 
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Proof. By Lemma -9ft(A(/3)) > for any f3 G [0, 1], so Ai > 0. Suppose 
Ai = 0. Because the eigenvalue A(/3) is a continuous function of the matrix 
f((3), where f((3) is also continuous on [0, 1] with respect to (3, so A(/3) is a 
continuous function of f3 on [0, 1]. This is due to the fact that a composition 
of continuous functions is still continuous. Noticing is also a continuous 
function, so — 9?(A(/3)) is continuous with respect to A(/3), and thus with 
respect to /3 on [0, 1]. Since a continuous function on a closed interval can 
achieve its minimum, there exists (3q G [0, 1] such that — 3R(A(/3o)) achieves 
the minimum Ai, i.e. — 3R(A(/3o)) = Ai = 0. This is contradicted to Lemma[71 
Therefore, Ai > 0. □ 

For any t G [0, oo), 0(t) G [0, 1]. Noticing B{t) = f(/3(t)), therefore 

{A | A is JB(t)'s non-zero eigenvalue, t > 0} 
={A | A is /(/3(t))'s non-zero eigenvalue, t > 0} (79) 
C{A | A is /(/9)'s non-zero eigenvalue; (3 G [0, 1]}. 



(80) 



Thus, 

A = inf{— 3ft(A) | A is _B(t)'s nonzero eigenvalue, t > 0} 

> inf{— 3?(A) | A is /(/3)'s non-zero eigenvalue; G [0, 1]} = A x . 

Because Ai > by Lemma [HJ so A > 0. That is, 
Corollary 1. Let 

A = inf{— 3?(A(t)) | X(t) is B(t) 's non-zero eigenvalue}, 

then A > 0. 

Based on this corollary, we can prove the approximation of X(t) to X(t). 

Proposition 12. Let X(t) = (xi(t), x 2 {t)) T = e tB{t) X(0), i.e. the solution 
of (E2jj. Let X(t) be defined by i.e. 

^■>=(W>)- u «( l o ) ( SIS ) ■ (81) 

Then lim \\X(t) - X(t) \\ = 0. 

t— >oo 
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Proof. Notice that eigenvectors of a matrix are continuous functions of the 
matrix. Since g(/3) is composed of the eigenvectors of the matrix f(0) and 
f(0) is continuous on [0, 1] with respect to 0, therefore g(0) is continuous 
on [0, 1] with respect to 0. Because the inverse of a matrix is a continuous 
mapping, so g~ 1 (0), i.e. the inverse of g(/3), is continuous with respect to 
g(0), and therefore is continuous on [0,1] with respect to since g(/3) is 
continuous on [0, 1]. Since any continuous function is bounded on a compact 
set [0,1], both g{0) and g~ 1 (0) are bounded on [0, 1]. That is, there exists 
K > such that \\g(/3)\\ < K and < K for a11 P £ [°, Because 

{U(t) \ te [0,oo)}= {g(J3(t)) | t G [0, oo)} C {g(0) | G [0, 1]}, 



we have 



sup ||*7(t)|| < sup 

*>o /?e[o,i] 



< K. 



Similarly, sup llf/ 1 (t)\\ < K. Notice 

t>o 



X(t) - X(t) =e tB{t) X(0) 
U(t) 



-U(t) 



1 

e 


e *A(i) 



[/(*) 





tA(t) 



c/(t)- 

t/(t)^X(0) 



-U(t) 



U(t) 



-1 



X(0) 



where X(t) is _B(t)'s nonzero eigenvalue. By a similar argument to X(0), X(t) 
is also real. Therefore, 

-X(t) = A(t)) = -3fr(A(t)) > A > 

or A(t) < — A < 0, where A is defined in Corollary [U Then 

) 



\X{t)-X 



U(t) 

mm 





^ 


(o 


e U(t) 1 



II^HIITOII 



<i^ 2 

<K 2 



|AT(0)||e <AW 
|X(0)||e-* A . 



Here we have used the norm property: 
have lim \\X(t) -X(t)\\ = 0. 

t— >oo 



\AB\\ < \\A\\\\B\ 



Since A > 0, we 
□ 
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2: An lower-bound estimation on population in local derivatives 

In the following, we will prove that there exists T, such that X\(t) > cM 
for any t > T. We first define a function 

h(B) = (h l (B),h 2 (B)) T = 9(P)( l n VW 



J y y > 1 \ x 2 (0) 

Clearly, we have X(t) = h{B{t)) and x x (t) = h^Bit)). Since B{t) e [0, 1] for 
all t, the following proposition can imply Xi(t) > cM . 

Proposition 13. There exists c > such that 

inf hAB) > cM. 
/3e[o,i] 

where M = xi(0) + £2(0), c is independent of M. 

Proof. Without loss of generality, we assume M — 1. We will show inf hAB) 

/3e[o,i] 

c > 0. Since is a continuous function of B which is due to the continuity 
of g(B) and g~ l (B), h\(B) can achieve its minimum on [0, 1]. That is, there 
exists Bq e [0, 1], such that 

/11(A)) = inf hx(P) = c. 
/9e[o,i] 



Consider the matrix 



and a set of linear ODEs 



di 
di 



— aBo b 
aB —b 



The solution of ( 152]) . given an initial value Z(0) = X(0) = (xi(0), x 2 (0)) T , is 
Z(t) = e^)X(O). 

According to (1711) . /(A)) can be diagonalised as 

HM = 9(M( ^ o) )9-'(A). 
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where A(/3 ) is the nonzero and real eigenvalue of f((3o). Thus 

Z(t) = e*'™X(p) = g%) ( I J [M ) s(A))- 1 ( ) • (83) 

Because A(/3 ) < by Lemma so as time goes to infinity, 



m = 9(M ( I J m ) Mr 1 ( 



(84) 



That is, lim Z(t) = h((3o). In the following, we discuss two possible cases: 

t— >oo 

I3 >0 and /3 = 0. 

If /?o > 0, then the transpose of the matrix f((3 ), i.e. f((3o) T , is an 
infinitesimal generator of an irreducible CTMC, which has two states and 
the transition rates between these two states are a/?o and b respectively. 
Moreover, the transient distribution of this CTMC, denoted by Z(t) = 
(zi(t), Z2(t)) T , satisfies the ODEs (18"2|) . As time goes to infinity, the tran- 
sient distribution Z(t) converges to the unique steady-state probability dis- 
tribution. Since lim Z(t) = h(/3o), therefore h((3o) = (hi(Po), /i2(A))) T is the 

t— >oo 

steady-state probability distribution and thus h\(Po) > 0. So inf^ 6 [ 0il ] h\{0) = 
MA>) > 0. ' 
If p = 0, then 

' b 







— — , — — I converges to zero. Letting 
at at J 

time go to infinity on the both sides of fl82|) . we obtain the following equilib- 
rium equations, 



oJ-^{w)- (85) 

By the conservation law, (foi(A))+^2(A))) T — M — 1. Therefore, (hi(Po), /i2(A))) J 

Sett lsflGS 

/(/3 )(/ il (/3o),^(/3o)) T = 0, 



h 1 (p Q ) + h 2 (Po) = M = l. (86) 

Solving fl86|) . we obtain the unique solution (hi(f3 ), h 2 ((3o)) T — (1, 0) T . There- 
fore, hi(Po) is one, and thus infg e [ 0i i] hi(/3) = hi(fto) > 0. □ 
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Remark 4. As (3 tends to 0, 



a/3 -b J Jy 1 \ -b 

Correspondingly, for the equilibrium h{/3) = (hi(/3), h 2 {f3)) T = (jj^j;, ^fff ) 

satisfying f{f3)h{f3) = and h x {p) + h 2 {(3) = M, we have ( J^, ggQ 
(M, 0) T as (3 tends to zero. From the explicit expression, i.e. h\(/3) - 



a/3+b' 



the minimum and maximum ofhi(/3) are and M respectively, which cor- 
respond to the matrices /(l) and /(O) respectively. In the context of the 
PEPA model, /(l) corresponds to a free subsystem and there is no synchro- 
nisation effect on it, i.e. the subsystem of component type X is independent 
of Y . The matrix /(O) reflects that the subsystem of X has been influenced 
by the subsystem ofY, i.e. the rates of shared activities are determined by Y , 
that is, the term a min{a:i, y±} has been replaced by ay\. Therefore the exit 
rates from the local derivative X± correspondingly become smaller since now 
ayi < ax\ . In order to balance the flux, which is described by the equilibrium 
equation, the population of X\ must increase. That is why the equilibrium 
hi(/3) increases as (3 decreases. In short, synchronisations can increase the 
populations in syncrhonised local derivatives in the steady state. 

As an application of the above facts, if hi((3o) > for some (3 > 0, then 
we can claim that hi(0) > because hi(0) > hi((3 ) > 0. 

Obviously, Proposition [13] has a corollary: 
Corollary 2. There exists c > such that for any t G [0, oo), x±(t) > cM. 

Proposition [T5] and Proposition [T3] can lead to the following lemma. 

Lemma 9. There exists c > 0, T > 0, such that x±(t) > cM for all t > T. 

Proof. By Proposition [TBI or Corollary El there exists ci,Xi > such that 
xi(t) > C\M for any t > T\. By Proposition [12], lim^oo \x\{€) — X\{€)\ = 0, 
which implies that for any e, there exists T 2 > such that for any t > T 2 , 
Xi(t) > Xi(t) — e. Choose T 2 > T 1; then we have 

2?i(t) > Xi(t)-e> dM - e. 

Therefore, there exist c, T > such that such that Xi(t) > cM for all t > 
T. □ 
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Because xx(t) > cM, provided cM > N we have xi(t) > cM > N > yi(t), 
i.e., the system will finally become linear. In the following we will show how 
to apply our method to more general cases. 

7.2. Convergence for two component types and one synchronisation (II): gen- 



This section deals with such an arbitrary PEPA model which has two 
component types and one synchronisation. The local action rates of the 
shared activity are not assumed to be the same. The main result of this 
section is a convergence theorem: as long as the population of one component 
type is sufficiently larger than the population of the other, then the solution 
of the derived ODEs converges as time tends to infinity. 

7.2.1. Features of coefficient matrix 

We assume the component types to be X and Y . The component type 
X is assumed to have local derivatives Xi,X2, • • ■ ,X m , while Y has local 
derivatives Yx,Y 2 , ■ ■ ■ ,Y n . We use Xi(t) to denote the population of X in 
Xj (i = 1, • ■ ■ ,m) at time t. Similarly, yj(t) denotes the population of Y 
in Yj (j = 1, ■ • • ,n) at time t. Without loss of generality, we assume the 
synchronisation is associated with the local derivatives X\ and Y\, i.e. the 
nonlinear term in the derived ODEs is min{rxi(t), syi(t)} where r and s are 
some constants. In fact, if the synchronisation is associated with Xj and Yj, 
by appropriately permuting their suffixes, i.e. i — > 1, i + 1 — > 2, ••• ,i — 1 — > 
m, j — > 1, j + 1 — >■ 2, • • • , j — 1 — >■ n, the synchronisation will be associated 
with X\ and Y\. According to the mapping semantics presented previously, 
the derived ODEs from this class of PEPA model are 



where x = (xi(t), ■ ■ ■ , x m (t), yi(t), • • • , y n {t)) T . For convenience, we denote 



In (1ST)) all terms are linear except for those containing "minjrxi (t) , syi (t) }" . 
Notice 

min{ra;i(t),sj/i(t)} = I{rxi{t)<a V1 (t)}rx 1 (t) + i{rai(t)>«i/i(t)}S2/i(*)- 



eral case 




(87) 



X(t) = {x l {t),x 2 {t),--- ,x m {t)) T , 
Y(t) = (yx(t),y 2 (t),--- ,y n (t)) T . 
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Whenrxi(t) < sy^t), which is indicated by I{rxi(t)<syi(t)} = 1 and I{ rxi (t)>syi(t)} 
0, we can replace min{rxi(t), syi(t)} by rxi(t) in (jSTJ). Then (I8TI) becomes 
linear since all nonlinear terms are replaced by linear terms rxi(t), so the 
ODEs have the following form, 



AX 

At 
AY 

At 



Y 



where Qi is a coefficient matrix. Similarly, if rx\{t) > syi(t), min{rxi(t), sy\(t)} 
can be replaced by syi(t) in ( 1H71) . Then ( 1H7I) can become 

f 1 = Qi ( y V ( 89 ) 



where Q2 is another coefficient matrix corresponding to the case of rx\(t) > 

syi(t). 

In short, the derived ODEs (IBT1) are just the following 



( it J = hrxi<s Vl }Ql ( y j + / {r-a:i> S2 /i}Q2 ( y J ■ (90) 

The case discussed in the previous section is a special case of this kind of 
form. If the conditions rxi(t) < syi(t) and rx\(t) > syi(t) occur alternately, 
then the matrices Qi and Q2 will correspondingly alternately dominate the 
system, as Figure [7] illustrates. 

Similar to the cases discussed in the previous two sections, the conver- 
gence problem of ( )90l) can be divided into two subproblems, i.e. to examine 
whether the following two properties hold: 

1. There exists a time T, such that either x\ < y±, Vt > T or x\ < 

Vt > T. 

2. The eigenvalues of Q\ and Q2 other than zeros have negative real parts. 

The first item can guarantee ( 190]) to eventually have a constant linear form, 
while the second item ensures the convergence of the bounded solution of 
the linear ODEs. If the answers to these two problems are both positive, 
then the convergence of the solution of (19"U|) will hold. The study of these 
two problems are discussed in the next two subsections. In the remainder of 
this subsection, we first investigate the structure property of the coefficient 
matrices Qi and Q2 in f l90|) . 
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(dX\ 




V dt J 



Figure 7: Illustration of derived ODEs with component types and one synchronisation 



The structure of the coefficient matrices Q± and Q2 is determined by the 
following two propositions, which indicate that they are either block lower- 
triangular or block upper-triangular. 



Proposition 14. Q\ in can be written as 

Qi = 



Qi 

w v nxn 



(91) 



(m+n)x (m+n) 



where Qj is an infinitesimal generator matrix with the dimension m x m, 
and 

( wu • • ■ \ 
w 2 i ••■ 



(92) 



\w nl ■ • • J 

where Wu < 0, Wji(j = 2, • • • ,n) > and Y^=i w ji = 0- Here V and W 
satisfy that if we let 



P = (W 1 + V ll V 2 . 



v n ) 



(93) 
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i. e. P 's first column is the sum of V 's first column and W 's first column, 
while P 's other columns are the same to V 's other columns, then P T is also 
an infinitesimal generator matrix. 



ft = ( t l\. 



Proof. Let 

a, = i 

v w v 

where Q\ and V are mxm and n x n matrices respectively. Suppose rx\{t) < 
syi(t), then 

y){y)- (95) 

d X 

— = Q l X + UY. (96) 

The condition rxi(t) < syi(t) implies that all nonlinear terms mm{rxi(t), syi(t)} 
can be replaced by X\(t). This means that the behaviour of the component 
type X in (1931) and (1931) is independent of the component type Y . Thus 
in (1931) [/ must be a zero matrix, i.e. 



dt 



So we have 



Qi 



Moreover, (1931) becomes 



Qi o 
w v 



IT = (97) 



that is, there is no synchronisation in the ODEs corresponding to the com- 
ponent type X given rxi(t) < syi(t). Then by Proposition HI is an 
infinitesimal generator. 
According to (193j) . 



dF 

_ = iyx + VY 
dt 



(W h W 2 ,--- ,W m )(x h x 2 ,--- ,x m y +VY 



(98) 



X^ + VY + ^XiWi.. 

i=2 
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where W = (Wi, W 2 , • • • , W m ). Notice that the component type Y is syn- 
chronised with the component type X only through the term mm{rxi(t), syi(t)} = 
x\(t). In other words, in fl98l) Y is directly dependent on only x\ other than 
Xi (i > 2). This implies W\ 7^ while Wi = (i — 2, 3, • ■ ■ , m). Therefore, 



where Vj (j = 1, 2, ■ • • , n) are the columns of V. Denote W\ = (wn, to 2 i, ■ ■ • , w n i) T . 
Notice that Y x is a pre local derivative of the shared activity, and Xiiu u rep- 
resents the exit rates of the shared activity from Y\. Therefore, W\\ < 0. 
Moreover, x±Wji (j — 2, • • • , n) are the synchronised entry rates for the local 
derivatives Yj (j = 2, ••• ,n) respectively, so Wji > (j = 2, • • • ,n). By 
the conservation law, the total synchronised exit rates are equal to the total 
synchronised entry rates, i.e. x\ YTj=i w ji = or YTj=\ w ji = 0- 

We have known that x\ in (|99|) derives from the synchronised term min{rxi, syi}. 
If the effect of the synchronisation on the behaviour of Y is removed, i.e. re- 
cover yi by replacing x\, then (|9"9"|) will become 



where P = (Wi + Vi, V2, ■ ■ ■ , V n ). Since there is no synchronisation contained 
in the subsystem of the component type Y, according to Proposition HI P T 




(99) 




(100) 



is the infinitesimal generator. 



□ 



Similarly, we can prove 



Proposition 15. Q2 in / fgQj) can be written as 




(101) 
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where fu < 0, fji(j = 2, ••• ,m) > and J2f=i fji = 0- Here F and E 
satisfy that if we let 

R = {F 1 + E 1 ,E 2 ,--- ,E n ), (103) 
then R T is also an infinitesimal generator matrix. 

7.2.2. Eigenvalues of coefficient matrix 

In this subsection, we will determine the eigenvalue property of Qi and 
Q2. First, the Perron- Frobenius theorem gives an estimation of eigenvalues 
for nonnegative matrices. 

Theorem 6. (Perron- Frobenius) . Let A = (a^) be a real n x n matrix 
with nonnegative entries > 0. Then the following statements hold: 

1. There is a real eigenvalue r of A such that any other eigenvalue A 
satisfies |A| < r. 

2. r satisfies min a^- < r < max £\ a «i- 

Remark 5. We should point out that in the second property, exchanging i 
and j in aij in the formula, we still have min J? , < r < max^ .a^. In 

fact, A T is also a real matrix with non-negative entries. Since A T and A 
share the same eigenvalues, so r is one of the eigenvalues of A T , such that 
any other eigenvalue A of A T satisfies |A| < r. Notice {A T )ij = Aj i; By 
applying the Perron- Frobenius theorem to A T , we have 

min dji < r < max a^. 

j j 

We cannot directly apply this theorem to our coefficient matrices Q\ and 
Q2, since both of them have negative elements, not only on the diagonal but 
also in other entries. However, we use some well-known techniques in linear 
algebra, i.e. the following Lemma [10] and [11] (which can be easily found 
in linear algebra textbooks), to cope with this problem, and thus derive 
estimates of their eigenvalues. 

Lemma 10. If E mxm and F nxn have eigenvalues Ai (i — 1, 2, ■ • ■ , m) and 5j 
(j = 1, 2, • ■ ■ , n) respectively, then each of 

Hi = (g f)'^ 2= (o f) 

has eigenvalues Aj (i — 1, 2, • • • , m) and 5j (j = 1, 2, • • ■ , n). 
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Lemma 11. If X is an eigenvalue ofV, then X + r is an eigenvalue ofV + rl, 
where r is a scalar. 



Theorem 7. The eigenvalues of both Q\ and Q2 are either zeros or have 
negative real parts. 

Proof. We only give the proof for Q±s case. By Proposition (HJ 

Qi=(i 1 , J ) • (104) 



w v n> 



(m+n) x (m+n) 



According to Lemma [TOj if all eigenvalues of Q\ and V are determined, then 
the eigenvalues of Q\ will be determined. Let us consider V first. 

Notice that only diagonal elements of V are possibly negative (which can 
be deduced from Proposition [T4"|) . Let r = sup \ Vu\ > 0, then all the entries 

i 

of V + rl are nonnegative. Let A be an arbitrary eigenvalue of V, then by 
Lemma [TT], A + r is an eigenvalue of V + rl. 

Notice the sum of the elements of any column of V is zero (because the 
sum of entry rates equals to the sum of exit rates) , so V has a zero eigenvalue 
with the corresponding eigenvector 1, i.e. VI = 0. Thus r = + r is an 
eigenvalue of V + rl. Moreover, 



in \ (V + rl)ji = r = max % (V + rl)^. 



mm 



Applying the Perron-Frobenius theorem (Theorem[H]) and RemarkOto V+rl, 
so 

\X + r\<r. (105) 

Let A = a + bi, then (I105P implies that a < 0, and if a = then b = 0. In 
other words, V^'s eigenvalues are either zeros or have negative real parts. 

Similarly, Q±s eigenvalues other than zeros have negative real parts. By 
Lemma [TO], the eigenvalues of Q\ other than zeros have negative real parts. 
The proof is complete. □ 

7.2.3. Convergence theorem 

Now we deal with another subproblem: whether or not after a long time, 
we always have rx\ > syi (or rx\ < sy\). If the population of X is signif- 
icantly larger than the population of Y, intuitively, there will finally be a 
greater number of X in the local derivative X\, than the number of Y in Y\. 
This will lead to rx\ > sy\. 
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Lemma 12. Under the assumptions in Section 7.1.1, for the following ODEs 
_ T n ( X\ T _ fx \ 

dK I _ ^{rxtKsy^Ql I y J + -*{ra;i>sj/i} V2 I y I , 

t/iere exists c\ > 0,C2 > 0,T > 0, stzc/i i/itrf £i(t) > ciM, j/i(t) > C2N for 
any t >T, where c\ and C2 are independent of M and N. 

Proof. The proof is essentially the same as the proof of Lemma [H We only 
give the sketch of the proof for X\(t) > C\M. By introducing new two func- 
tions a{t) and (3(t), 

( min{r gl (tW(t)} x At) ^ 0, 

1, xi(t)=0, 



1 



= - a(s)ds, 
t Jo 

the nonlinear term min{rxi(i), s?/i(t)} equals ra(t)xi(t), and thus the ODEs 
associated with the subsystem X can be written as 

H X 

— = A(t)X (106) 

where A{t) is related to a{t). The solution of fTT06ll is X(t) = e tB(t) X(0), 

1 /•* 

where B(t) is defined by B(t) — - I A(s)ds, and thus B(t) is related to 

t Jo 

0(t)- 

Notice that according to Theorem [7] and its proof, the eigenvalues of B(t) 
other than zeros have negative real parts for any t > 0. By a similar proof 
to Corollary [TJ we have 

A = inf{— 3?(A) | A is f?(t)'s non-zero eigenvalue} > 0. (107) 



i>0 



This fact will lead to the conclusion that X(t) can be approximated by X(t), 
where X(t) is constructed similarly to the one in Proposition [12j Because a 
general B(t) considered here may not be diagonalisable, so the construction 
of X(t) is a little bit more complicated. We detail the construction as well 



as the proof of the following result in Appendix D 



lim \\X(t) - X(t)\\ = 0. 

t— >oo 
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Then, by similar arguments to Proposition [T3] and Corollary 121 we can prove 
that inf t> jn Xi(t) > cM, where Xi(t) is the first entry of X(t). Then, by a 
similar proof to the proof of Lemma [9], we can conclude that there exists a 
number c\ such that x\(t) > c\M after a time T. □ 

Lemma 13. Under the assumptions of Lemma [7D[ if M > K\N or N > 
K 2 M, where constants K\ > and K 2 > are sufficiently large, then there 
exists T such that rx\{t) > syi(t), Vt > T or rx\{t) < syx(t), Wt > T 
respectively. 

Proof. By the boundedness of solutions, < Xi(t) < M and < y\{t) < N 
for any t. Suppose M > KiN . By Lemma [T2l there exists c, T > 0, Xi(t) > 
cM, Vt > T. Since K\ is assumed to be large enough such that K\ > 
then rcM > sN. So for any t > T, we have 

rxiit) > rcM > sN > sy^t). 

If N > K 2 M, the proof is similar and omitted here. □ 

Now we state our convergence theorem. 

Theorem 8. If M > K\N or N > K 2 M, where constants Ki, K 2 > are 
sufficiently large, then the solution of the derived ODEs ( fPOj) . i. e. 

^ \ _ T _ / X\ ( X\ 

J — Hrx^syxjQl I y J Hrx 1 >sy 1 }Q2 I y I , 

converges to a finite limit as time goes to infinity. 

Proof. Suppose M > K\N, then by Lemma [131, there exists a time T > 0, 
such that after time T, rx\{t) > syi(t), so ( 1901) becomes 

f ) = Q* ( y V (108) 



dt 



Since (52's eigenvalues other than zeros have strict negative real parts ac- 
cording to Theorem [71 and the solution of the above equation is bounded, 
then by Corollary [31 the solution converges to a finite limit as time goes to 
infinity. □ 
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Remark 6. We should point out that for a general PEPA model with two 
component types and one synchronisation, the limit of the solution of the 
derived ODEs is determined by the populations of these two component types 
rather than the particular starting state. That is to say, whatever the initial 
state is, as long as the total populations of the components are the same, the 
limit that the solution converges to will always be the same. We do not plan 
to discuss this topic in detail in this paper. 



8. Related work 



As we have mentioned in the introduction section, the fluid approxima- 
tion approach was first proposed by Hillston in 2l| to deal with the state 
space explosion problem in the context of PEPA. An interpretation as well 
as a justification of this approximation a ppr oach has been demonstrated by 
Hayden in his dissertation 



33 . In 33 



|, generation of similar systems 
of coupled ODEs for higher-order moments such as variance has been ad- 
dressed. Additionally, the dissertation 33| discusses how to derive stochastic 
differential equations from PEPA models. 

More recently, some extensions of the previous mapping from PEPA to 
ODEs have been presented by Bradley et al. in [23|. In particular, passive 



rates are introduced into the fluid approximation. In the recent paper [35 



different existing styles of passive cooperation in fluid models are compared 
and intensively discussed. Moreover, a new passive fluid semantics for passive 
cooperation, which can be viewed as approximating the first moments of the 
component counting processes, has been provided, with a theoretical justifi- 



cation. The paper 36] considers the application of this fluid approximation 
approach with modifications in the context of epidemiology. In this paper, 
the notions of side and self-loops are added to the activity matrix, and the 
rates are calculated differently, for the purpose of deriving from PEPA mod 
els the most commonly used ODEs in the context of epidemiology. In 
and 
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38 



by Tribastone, a new operational semantics is proposed to give a 
compact symbolic representation of PEPA models. This semantics extends 
the application scope of the fluid approximation of PEPA by incorporating 
all the operators of the language and removing earlier assumptions on the 
syntactical structure of the models amenable to this analysis. 

The fluid approximation approach has also been applied to timed Petri 



nets to deal with the state space explosion problem [39|,|40|. The comparison 



between the fluid approximation of PEPA models and timed continuous Petri 
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nets has been demonstrated by Galpin in 41]. This paper has established 



links between two continuous approaches to modelling the performance of 
systems. In the paper, a translation from PEPA models to continuous Petri 
nets and vice versa has been presented. In addition, it has been shown 
that the continuous approximation using PEPA has infinite server semantics. 
The fluid approximation approach has also been used by Thomas to derive 
asymptotic solutions for a class of closed queueing networks 42j. In this 
paper, an analytical solution to a class of models, specified using PEPA, 
is derived through the ODE approach. It is shown that "this solution is 
identical to that used for many years as an asymptotic solution to the mean 
value analysis of closed queueing networks" . 

Moreover, the relationship between the fluid approximation and the un- 
derlying CTMCs for a special PEPA model has been revealed by Geisweiller 
et al. in 25]: the ODEs derived from the PEPA description are the limits 
of the sequence of underlying CTMCs. It has been shown in 28j by Gilmore 
that for some special examples the equilibrium points of the ODEs derived 
from PEPA models coincide the steady-state probability distributions of the 
CTMCs underlying the nonsynchronised PEPA models. 

In addition, there are several papers which discuss how to derive response 



time from the fluid approximation of PEPA models. In (43|, by constructing 
an absorption operator for the PEPA language, Bradley et al. allow gen- 
eral PEPA models to be analysed for fluid-generated response times. Clark 
et al. demonstrate in 44( how to derive expected passage response times 
using Little's law based on averaged populations of entities in an equilib- 
rium state. This technique has been generalised into one for obtaining a full 
response-time profile computing the probability of a component observing the 
completion of a response at a given time after the initiation of the request, 



sec 
in 



45 



Moreover, an error in the passage specification in the approach taken 
43| has been uncovered and rectified in 45 1 by Clark. The ODE method 



associated with the PEPA language has demonstrated successful application 



in the performance analysis of large scale systems [46|, |47|, |48|, |49|, |50 



9. Conclusions 

In this paper, we have demonstrated how to derive the fluid approxi- 
mation from a general PEPA model via the numerical representation of the 
PEPA language, which extended the current mapping semantics of fluid ap- 
proximations. The fundamental properties of the fluid approximation such 
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as the existence and uniqueness, boundedness and nonnegativeness of the 
solutions of the derived ODEs have been established. Moreover, the con- 
vergence of the solution under a particular condition for general models has 
been verified. This particular condition relates some famous constants of 
Markov chains such as the spectral gap and the Log-Sobolev constant. For 
the models without synchronisations or with one synchronisation and two 
component types, we have determined the convergence. Furthermore, the 
consistency between the fluid approximation and the CTMC in the context 
of PEPA has been revealed. If a model has no synchronisations, then the 
derived ODEs are just the probability distribution evolution equations of the 
underlying CTMC except for a scaling factor. For any general PEPA model, 
the ODEs can be taken as the corresponding density dependent CTMC with 
the concentration level infinity. In addition, the coefficient matrices of some 
derived ODEs were studied: their eigenvalues are either zeros or have nega- 
tive real parts. The structural property of invariance has been shown to play 
an important role in the proof of convergence for some PEPA models. Due 
to the limitation of pages, we have not shown how to derive performance 
measurers from the fluid approximation, and have not given the numerical 
comparison between the fluid approximation and the CTMC in terms of 



performance measures. For more details, please refer to [22 



In addition to the established results, this paper has also demonstrated 
comprehensive techniques and methods to investigate the PEPA language: 
not only both theoretical and experimental, probabilistic and analytic, but 
also syntactic and numerical (investigating the models based on the numerical 
representation of PEPA), qualitative and quantitative (exploiting the struc- 
tural property to verify the convergence). These results and techniques are 
expected to have more applications in performance modelling and evaluation 
of large scale systems. 
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Appendix A. 



Theorem 9. (Kurtz theorem 13W . page 456). Let X n be a family of density 
dependent CTMCs with the infinitesimal generators 

9*3+1 = nf(k/n,l), 

where f(x, I) (x G E C R h , I G Z' 1 ) is a continuous function, k is a numerical 
state vector and I is a transition vector. 
Suppose X(t) G E satisfies 

dx „. , 



where F(x) = ^ //(x, /). Suppose that for each compact K C J5 ; 

V" ||/|| Bup/(z,0 < oo 

and t/iere exzsis > snc/i £/iat 

||F(x)-F(y)|| <M K \\x-y\\, x,yeK. 
= x , then for every t > 0, 



lim sup 



s<t 



n 



-X(s) 



a.s. 



(A.l) 



(A.2) 



(A.S) 



The following lemma can be found in any good book on differential cal- 
culus. 

Lemma 14. Let y(t) be a differentiate function defined for t > 0. Suppose 

dy 

a, b G R, d^O. If y(t) satisfies — > ay(t) + 6, t > 0, then 

Lit- 



y(t) > e at ( y(0) + -)--. 

a / a 



dy 

Similarly, if y(t) satisfies — < ay(t) + b, t > 0, then 



y(t) < e at [ y(0) + -)--. 

a / a 
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Proof. Let W{t) = y(t)e~ a \ then 



^ = e~ at - ay 



ing on both sides, so W(t)-W(fS) > b J e" as ds. Thus y(t)e~ at -y(0) > J(l- 
e~ at ). 

So > e a * (2/(0) + |) — ^. The second conclusion can be similarly proved. 

□ 



Theorem 10. (Fundamental Inequality, f^jj . page 14)- If —r- = f(x, t) 

is defined on a set U inM. n x R with the Lipschitz condition 

||f( Xl ,t)-f(x 2 ,t)|| < J fr||x 1 -x 2 || 

for all (xi,t) and (x 2 ,t) on U , and if for €i,5 £ R 7 and Ui(i) and u 2 (i) are 
£u>o continuous, piecewise differentiate functions on U into W 1 with 



dui(t) 



df 



f(u,(t),t) 



< e*, ana" ||ui(t ) - u 2 (t )|| < 5, 



then 



\Mt) - u 2 (t)|| < + {^p^j {e K{t - to) ~ 1) 



Appendix B. 



The material presented here is extracted from 32|]. Let (K,ir) be a 
Markov chain on a finite set S, where K is a Markov kernel and it is the 
stationary probability distribution associated with K. For any real function 
/, g on S, define an inner product "(•, •)" as 



(/,#} = ^2f(x)g(x)ir(x). 



Denote 



v / (/>/) ) and 

1*00 = {/ = 



< 00}. 



Then l 2 (n) is a Hilbert space with the norm || ■ || 2 . We say that K* is adjoint 
to K if 

(Kf,g) = {f,K*g), Vf,g £ L 2 (vr). 
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It follows that 

K*(x,y) = ^\K(y,x). 
ir[x) 

\i K = K*, then K is called self- adjoint. If K is self-adjoint on I 2 (it), then 
(K, ir) is reversible. 

For a function in / G l 2 {^), denote its mean and variance by 7r(/) and 
Var(f) respectively, that is 

= E ^ ar (/) = *((/ - *(/)) 3 )- 

Definition 8. (Dirichlet form). The form 

S(f,g) = ((I-K)f,g} 
is called the Dirichlet form associated with H t = e - *^ - ^ . 
Remark 7. T7ie Dirichlet form £ satisfies 

£(f, f) = ((/ - K)f, f) = ((l- ^^-) /, / 



s{fj) = \Y J U^)-f{y)) 2 K{x^{x). 



X.IJ 



Definition 9. (Spectral gap). Let K be a Markov kernel with Dirichlet 
form £ . The spectral gap A = X(K) is defined by 

Remark 8. In general A is the smallest non zero eigenvalue of I — K+ 2 K * ■ 
If K is self-adjoint, then A is the smallest non zero eigenvalue of I — K . 
Clearly, we also have 

A = min{£(/,/):||/|| 2 = l,7r(/) = 0}. 

The definition of the logarithmic Sobolev (Log-Sobolev) constant a is 
similar to that of the spectral gap A where the variance has been replaced by 



£(/) = £/(*) a iog (jj£) < x )- 

XdzS 
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Definition 10. (Log-Sobolev constant). Let K be an irreducible Markov 
chain with stationary measure tt. The logarithmic Sobolev constant a = a(K) 
is defined by 

a=min {^7r 

Lemma 15. For any finite Markov chain K with stationary measure tt, the 
Log-Sobolev constant a and the spectral gap A satisfy 

1-2-kM ^ A 
-A < a < 



log[l/7r(*) - 1] " " 2' 
where tt(*) = min x n(x). 



Appendix C. 

In this appendix, we claim that all eigenvalues of Qi (i = 1,2,3,4) ap- 
pearing in ( 154")) in Section 16.2.1 1 other than zeros have negative real parts. We 
do not worry about Q\ and Q4, since they are lower or upper block triangular 
matrices and the eigenvalues of this kind of matrices can be well estimated: 
all eigenvalues of Q\ and Q4 are either zeros or have negative real parts. 
All that we want to do here is to show that both Q2 and Q3 also have this 
property. 

By symbolic calculation using Matlab, Qs's eigenvalues are 
Ai,2,3 = 0(three folds), A4 = — C4 — C3, 

A 5 = + °2 + c i + °2) + 7^V( a i ~ a 2 + c i + c 2 ) 2 - 4aic 2 , 

^6 = ~^( a i + a 2 + Ci + c 2 ) - - \/(ai - a 2 + Ci + c 2 ) 2 - 4aic 2 . 

If (ai — a 2 + C\ + c 2 ) 2 — 4a! c 2 < 0, then the real parts of A5 and \$ are 
~\{ a i + a 2 + c i + c 2 ), which is negative. Otherwise, 

(ai - a 2 + ci + c 2 ) 2 - 4a x c 2 > 0. 

In this case, 

(ai - a 2 + ci + c 2 ) 2 - 4aic 2 < (a x - a 2 + ci + c 2 ) 2 < (ai + a 2 + c 1 + c 2 ) 2 , 



SO 



so 



--(ai + a 2 + ci + c 2 ) + - - a 2 + ci + c 2 ) 2 - 4aic 2 < 0. 

This means that A5 and A6 are both negative. Thus A, (i = 1, 2, • • • , 6) are 
either or have negative real parts. 

Similarly, Q 2 's eigenvalues are Si 23 = 0, £4 = — c x — c 2 , 



£5 = -77(0-1 + a 2 + c 3 + c 4 ) + ^A/( fl 2 - «i + c 3 + C4) 2 - 4a 



2C3, 



5 6 = ~2^ ai + a 2 + c 3 + C4) - - ^/(^ - ai + c 3 + c 4 ) 2 - 4a 2 c 3 . 

By similar argument, we still have that 8i (i — 1, 2, • • • , 6) are either zeros or 
have negative real parts. 



Appendix D. 

Appendix D.l. 

In this subsection, we use 9ft(z) and to respectively represent the 
real and imaginary parts of a complex number z. The following is mainly 
extracted from 51[ (page 39~42). 

Theorem 11. (The Jordan Canonical Form). Let A be a real matrix 
with real eigenvalues Xj, j = 1, • • • , k and complex eigenvalues Xj = aj + ibj 



A, 



ibj, 



and 

j = k+1, ■ ■ ■ ,n. Then there exists a basis {vi, • • • , Vfc, v^+i, u^+i, • • • , v n , u n } 
for M? n ~ k where Vj, j = 1, • • • , k and Wj, j = k + 1, • • • , n are generalized 
eigenvectors of A, 

Uj = 3?(wj) and Vj = S(wj) for j = k + 1, • • • ,n, such that the matrix 
P = {vi, • • • , v fc , Vfc+i, Ufc+i, • • • , v n , u n } is mvertible and 



P- X AP 



B l 



B, 



where the elementary Jordan blocks B — Bj, j — 
form 

X 1 ••■ 



B 



A 

••• 
••• 



... 

A 1 
A 



r are either of the 



(D.l) 
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for A one of the real eigenvalues of A or of the form 



with 



D 





D 


h 





... o 







D 


h 


... o 


B -- 














• 






D I 2 




■ 






D 


a —b 




' 1 


" 


and 




, h = 







b a 




1 





(D.2) 








for A = a + ib one of the complex eigenvalues of A. 



The Jordan canonical form of A yields some explicit information about 
the form of x = e At x , i.e. the solution of the initial value problem 



dx 

— =Ax 
dt 



That is, 



x(0) =x 
x(t) = Pdiag [e B ^] P -1 x , 



where Bj are the elementary Jordan blocks of A, j 
represents 

••• \ 

o B 2t ... Q 

diag [e Bjt ] 



(D.3) 
(D.4) 

Here diag [e Sj *] 



/ e 



Bit 



\ 












-,B r t 



In the following, the notation diag[-] indicates the similar meaning. If Bj = B 
is an m x m matrix of the form (ID. ip and A is a real eigenvalue of A then 



1 









t t 2 /2\ 
1 t 
1 



t 



m— 1 



/{m 







1 





t m "7(m-2)! 
t m - 3 /(m-3)! 

t 
1 



(D.5) 
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If Bj = B is an 2m x 2m matrix of the form (1D.2|) and A = a + ib is a complex 
eigenvalue of A, then 



= e 



R Rt Rt 2 /2\ ■■■ Rt m - l /{m- 1)! 
i2 itt ••• Rt m ~ 2 /(m-2)\ 
i? ••• Rt m - 3 /(m-3)\ 



(D.6) 












Rt 
R 



where R is the rotation matrix 

^ _ cos bt — sin bt 
sin 6£ cos 6£ 

Theorem 12. 7/x(£) satisfies AD. Sty , then each coordinate in x(t) is a linear 
combination of functions of the form 



where X = a + ib is an eigenvalue of the matrix A nxn and < k < n — 1. 

Corollary 3. // £/ie eigenvalues of A are either zeros or have negative real 
parts, and x(£) is bounded in [0, oo), then x(i) converges to a finite limit as 
time goes to infinity. 

Proof. The solution is composed of the terms like t k e at cos bt and t k e at sin bt. 
If a < 0, then t k e at cosbt and t k e at sinbt converge as time goes to infinity. If 
a = b = 0, we will see k = 0. In fact, in this case t k e at cosbt = t k . If k > 0, 
then this term t k in the solution will make the solution unbounded as t tends 
to infinity. So k must be zero in the terms corresponding to a = b = 0. Thus, 
t k e at cosbt = 1 and t k e at sinbt = 0. So the solution converges. □ 

Appendix P. 2. 

As the above subsections illustrate, the following problem 



t k e at cosbt or t k e at sinbt 



( x(0) =x 
has solution x(£) = e A *x , which equals 

x(t) = Pdiag [e^*] P _1 xo, 




(D.7) 



S3 



where Bj are the elementary Jordan blocks of A, j = 1, • • • , r. Suppose the 
rank of A nxn is n — 1. This means that zero is a one fold eigenvalue of A. 
According to (1D.7|) . we construct a corresponding x in the form 



x = Pdiag Bj P _1 x , (D.8) 



where Bj(t) is defined as follows. If e Bjt in (1D.7j) has the form of f ]D.6j) . then 
Bj is defined by 

Bj = 2mX 2m- (D-9) 

If e Bjt has the form of (1D.5j) . and the corresponding real eigenvalue A < 0, 



then Bj is defined by 

Bj = O mxm (D.10) 

If A = 0, we know that zero is a one fold eigenvalue of A due to its rank 
n — 1. So m = 1. Then, 

4 = 1. (D.ll) 

In short, only for the zero eigenvalue is -Bj set to one, otherwise it is set to 
zeros. Clearly, we have 

Lemma 16. If zero is a one fold eigenvalue of A and all other eigenvalues 
of A have negative real parts, then 

lim |x(t) - x(t)| = 0. 

t— >oo 

Proof. If B is the Jordan block corresponding to the one fold zero eigenvalue, 
then according to (ID. lip . e Bt = e° = 1, and e Bt — 5=1 — 1 = 0. For any 
non-zero eigenvalue, since B — then e Bt — B = e Bt . Notice that by flD.5j) 
and (jEED , 

||e Bt || < C^e-^, 
where C\(t) is a polynomial of t with the maximum order k, and 

A = inf{— 3?(A) | A is non-zero eigenvalue of A} > 0. 

Therefore, 

||x(t) - x|| < W eBt ~ B W ^ C 2(t)e~ At — > (D.12) 
as t goes to infinity, where C^t) is a polynomial function of t. □ 
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The above construction can be extended to the problem of 

^§=A(t)X (D.13) 
with the initial value X(0), which is discussed in Section [7.2.31 The solution 



is 



X(t) = e tc{t) X(0), 



1 



where C(t) = — / A(s)ds. We can similarly define a function / such that 
t Jo 

C(t) = f((3(t)), where (3(t) is similarly defined according to A(t). Therefore, 

X(t) = e tc ^X(0) = e tf ^X(0). 

For a fixed /3, 

e tm X(0) = P(f3) diag [e B ^] P(/3)" 1 X(0), (D.14) 

where Bj(/3) are the elementary Jordan blocks of /(/?), j = 1, • • • ,r(f3). 
Repeating the previous construction process with Bj(j3) for each j, we obtain 
the constructed matrix B(f3)j. We define 

h(J3) = P(J3) diag \B(J3) 3 ] Pifi-'XiO). (D.15) 



For convenience, suppose the dimension of A(t) in (1D.13|) is n x n. We 



should point out that for any t, the rank of A(t) is n — 1, and thus for 
any t, A(t)'s zero eigenvalue is one fold. In fact, the rank of any infinitesimal 
generator with dimension nxn is n—1. This implies that any n— 1 columns or 
rows of this generator are linearly independent. According to the definition of 
A(t) in Section 17.2.31 A(t) is an infinitesimal generator if a(t) ^ 0. If a(t) = 
0, A(t) is also a generator after one column is modified (see Proposition 
and [15]), which means that the other n — 1 columns are linearly independent. 
So whatever t is, the rank of A(t) is n — 1. Thus, the zero eigenvalue is one 
fold. Therefore, /(/3)'s zero eigenvalue is also one fold for any /3. So each 
entry of all blocks B(/3)j is zero, except for the one corresponding to the zero 
eigenvalue, in which case this block is a scalar one. This implies that for any 
f3 all entries of the matrix diag [B (/3)j\ are zeros, except for a diagonal entry 
with one. 
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By permutation, diag B(/3)j can always be transformed into the form 



1 




Correspondingly, P(j3) is permuted into U((3). Therefore, the formulae (1D.15|) 
can be written as 



HP) = u(p) 



1 




(D.16) 



Now we prove a proposition which is used in Section 17.2.31 
Proposition 16. For the e*-^X(0) in ( |X>. 1$ and HP) in IIP. 15)) , we have 

lim \\e tm X(0)-h(/3)\\ =0. 

t— >co 

Proof. By a similar estimation as (ID. 121) . we have 



i tm X(0) -h{P)\\ < C(t)t 



-Ait 



where C(t) is a polynomial of t. By a similar proof to Lemma El we have 



Ai = inf ^{ — 3?(A)|A is /(/3)'s non-zero eigenvalue} > 0. 



Then 



MP)X(0) -h{P)\\ < C(t)e~ Alt 



as t goes to infinity. 



□ 



Let X(t) = h(P(t)), where P(t) G [0,1], and notice X(t) = f(P(t)). As a 
consequence of this proposition, we have 



Corollary 4. Let X(t) be the solution of ^t- = A(t)X which is discussed in 
Section YT2^ and let X(t) = h(P(t)), then 



lim \\X(t)-X(t)\\ = 0. 

t— >oo 
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Appendix D.3. 

This subsection presents a proof of Lemma [7] in Section 17.1.31 Let \(/3) 
be a nonzero eigenvalue of the following f(/3): 



We will prove the following lemma which states that the real part of A(/3) is 
negative. The proof given here does not rely on the explicit expression of the 
eigenvalue. 

Lemma [7J For any (3 e [0,1], 5J(A(/3)) < 0, where A(/3) is a nonzero 
eigenvalue of f(/3). 

Proof. After a shift max{a/3, b}I, f((3) becomes f(fl) = f(f3) + max{a/3, b}I, 
which is a nonnegative matrix. Then similarly to the proof of Theorem [71 
which is based on the Perron- Frobenious theorem (Theorem [6]), we can con- 
clude that the eigenvalue other than zero has negative real part. □ 
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